Line data Source code
1 : //# Copyright (C) 1997-2010
2 : //# Associated Universities, Inc. Washington DC, USA.
3 : //#
4 : //# This library is free software; you can redistribute it and/or modify it
5 : //# under the terms of the GNU Library General Public License as published by
6 : //# the Free Software Foundation; either version 2 of the License, or (at your
7 : //# option) any later version.
8 : //#
9 : //# This library is distributed in the hope that it will be useful, but WITHOUT
10 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
12 : //# License for more details.
13 : //#
14 : //# You should have received a copy of the GNU Library General Public License
15 : //# along with this library; if not, write to the Free Software Foundation,
16 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
17 : //#
18 : //# Correspondence concerning AIPS++ should be addressed as follows:
19 : //# Internet email: casa-feedback@nrao.edu.
20 : //# Postal address: AIPS++ Project Office
21 : //# National Radio Astronomy Observatory
22 : //# 520 Edgemont Road
23 : //# Charlottesville, VA 22903-2475 USA
24 : //#
25 : //# $Id: $
26 :
27 : #include <casacore/casa/Arrays/Matrix.h>
28 : #include <casacore/casa/Arrays/Cube.h>
29 : #include <casacore/casa/Arrays/ArrayMath.h>
30 : #include <casacore/casa/Arrays/MatrixMath.h>
31 : #include <casacore/casa/IO/ArrayIO.h>
32 : #include <casacore/casa/BasicMath/Math.h>
33 : #include <casacore/casa/BasicSL/Complex.h>
34 : #include <casacore/casa/Logging/LogIO.h>
35 : #include <casacore/casa/OS/File.h>
36 : #include <casacore/casa/Containers/Record.h>
37 :
38 : #include <casacore/lattices/LRegions/LCBox.h>
39 : #include <casacore/casa/Arrays/Slicer.h>
40 : #include <casacore/scimath/Mathematics/FFTServer.h>
41 : #include <casacore/casa/OS/HostInfo.h>
42 : #include <casacore/casa/Arrays/ArrayError.h>
43 : #include <casacore/casa/Arrays/ArrayIter.h>
44 : #include <casacore/casa/Arrays/VectorIter.h>
45 :
46 :
47 : #include <casacore/casa/Utilities/GenSort.h>
48 : #include <casacore/casa/BasicSL/String.h>
49 : #include <casacore/casa/Utilities/Assert.h>
50 : #include <casacore/casa/Utilities/Fallible.h>
51 :
52 : #include <casacore/casa/BasicSL/Constants.h>
53 :
54 : #include <casacore/casa/Logging/LogSink.h>
55 : #include <casacore/casa/Logging/LogMessage.h>
56 :
57 : #include <synthesis/MeasurementEquations/MatrixCleaner.h>
58 : #include <casacore/coordinates/Coordinates/TabularCoordinate.h>
59 : #ifdef _OPENMP
60 : #include <omp.h>
61 : #endif
62 : using namespace casacore;
63 : namespace casa { //# NAMESPACE CASA - BEGIN
64 :
65 :
66 0 : Bool MatrixCleaner::validatePsf(const Matrix<Float> & psf)
67 : {
68 0 : LogIO os(LogOrigin("MatrixCleaner", "validatePsf()", WHERE));
69 :
70 : // Find the peak of the raw Psf
71 0 : AlwaysAssert(psf.shape().product() != 0, AipsError);
72 0 : Float maxPsf=0;
73 0 : itsPositionPeakPsf=IPosition(psf.shape().nelements(), 0);
74 : //findMaxAbs(psf, maxPsf, itsPositionPeakPsf);
75 0 : Int psfSupport = findBeamPatch(0.0, psf.shape()(0), psf.shape()(1),
76 0 : 4.0, 20.0);
77 0 : findPSFMaxAbs(psf, maxPsf, itsPositionPeakPsf, psfSupport);
78 0 : os << "Peak of PSF = " << maxPsf << " at " << itsPositionPeakPsf
79 0 : << LogIO::POST;
80 0 : return true;
81 0 : }
82 :
83 :
84 0 : MatrixCleaner::MatrixCleaner():
85 0 : itsMask( ),
86 0 : itsSmallScaleBias(0.0),
87 0 : itsMaskThreshold(0.9),
88 0 : itsDirty( ),
89 0 : itsXfr( ),
90 0 : itsScaleSizes(0),
91 0 : itsMaximumResidual(0.0),
92 0 : itsStrengthOptimum(0.0),
93 0 : itsTotalFlux(0.0),
94 0 : itsChoose(true),
95 0 : itsDoSpeedup(false),
96 0 : itsIgnoreCenterBox(false),
97 0 : itsStopAtLargeScaleNegative(false),
98 0 : itsStopPointMode(-1),
99 0 : itsDidStopPointMode(false),
100 0 : itsJustStarting(true),
101 0 : psfShape_p(0),
102 0 : noClean_p(false)
103 : {
104 0 : itsMemoryMB=Double(HostInfo::memoryTotal()/1024)/16.0;
105 0 : itsScales.resize(0);
106 0 : itsScaleXfrs.resize(0);
107 0 : itsDirtyConvScales.resize(0);
108 0 : itsPsfConvScales.resize(0);
109 0 : itsScaleMasks.resize(0);
110 0 : itsScalesValid = false;
111 0 : itsStartingIter = 0;
112 0 : itsFracThreshold=Quantity(0.0, "%");
113 0 : }
114 :
115 :
116 :
117 0 : MatrixCleaner::MatrixCleaner(const Matrix<Float> & psf,
118 0 : const Matrix<Float> &dirty):
119 0 : itsMask( ),
120 0 : itsSmallScaleBias(0.0),
121 0 : itsScaleSizes(0),
122 0 : itsMaximumResidual(0.0),
123 0 : itsStrengthOptimum(0.),
124 0 : itsTotalFlux(0.0),
125 0 : itsChoose(true),
126 0 : itsDoSpeedup(false),
127 0 : itsIgnoreCenterBox(false),
128 0 : itsStopAtLargeScaleNegative(false),
129 0 : itsStopPointMode(-1),
130 0 : itsDidStopPointMode(false),
131 0 : itsJustStarting(true),
132 0 : noClean_p(false)
133 : {
134 0 : AlwaysAssert(validatePsf(psf), AipsError);
135 0 : psfShape_p.resize(0, false);
136 0 : psfShape_p=psf.shape();
137 : // Check that everything is the same dimension and that none of the
138 : // dimensions is zero length.
139 0 : AlwaysAssert(psf.shape().nelements() == dirty.shape().nelements(),
140 : AipsError);
141 0 : AlwaysAssert(dirty.shape().product() != 0, AipsError);
142 : // looks OK so make the convolver
143 :
144 : // We need to guess the memory use. For the moment, we'll assume
145 : // that about 4 scales will be used, giving about 32 TempLattices
146 : // in all. Also we'll try not to take more that half of the memory
147 :
148 : // Ah, but when we are doing a mosaic, its actually worse than this!
149 : // So, we pass it in
150 0 : itsMemoryMB=Double(HostInfo::memoryTotal()/1024)/16.0;
151 :
152 0 : itsDirty = new Matrix<Float>(dirty.shape());
153 0 : itsDirty->assign(dirty);
154 0 : setPsf(psf);
155 0 : itsScales.resize(0);
156 0 : itsScaleXfrs.resize(0);
157 0 : itsDirtyConvScales.resize(0);
158 0 : itsPsfConvScales.resize(0);
159 0 : itsScaleMasks.resize(0);
160 0 : itsScalesValid = false;
161 0 : itsStartingIter = 0;
162 0 : itsFracThreshold=Quantity(0.0, "%");
163 0 : }
164 :
165 :
166 0 : void MatrixCleaner::setPsf(const Matrix<Float>& psf){
167 0 : itsXfr=new Matrix<Complex>();
168 0 : AlwaysAssert(validatePsf(psf), AipsError);
169 0 : psfShape_p.resize(0, false);
170 0 : psfShape_p=psf.shape();
171 0 : FFTServer<Float,Complex> fft(psf.shape());
172 0 : fft.fft0(*itsXfr, psf);
173 : //cout << "shapes " << itsXfr->shape() << " psf " << psf.shape() << endl;
174 0 : }
175 :
176 0 : MatrixCleaner::MatrixCleaner(const MatrixCleaner & other)
177 :
178 : {
179 0 : operator=(other);
180 0 : }
181 :
182 0 : MatrixCleaner & MatrixCleaner::operator=(const MatrixCleaner & other) {
183 0 : if (this != &other) {
184 0 : itsCleanType = other.itsCleanType;
185 0 : itsXfr = other.itsXfr;
186 0 : itsMask = other.itsMask;
187 0 : itsDirty = other.itsDirty;
188 0 : itsScales = other.itsScales;
189 0 : itsScaleXfrs = other.itsScaleXfrs;
190 0 : itsPsfConvScales = other.itsPsfConvScales;
191 0 : itsDirtyConvScales = other.itsDirtyConvScales;
192 0 : itsScaleMasks = other.itsScaleMasks;
193 0 : itsStartingIter = other.itsStartingIter;
194 0 : itsMaximumResidual = other.itsMaximumResidual;
195 0 : itsIgnoreCenterBox = other.itsIgnoreCenterBox;
196 0 : itsSmallScaleBias = other.itsSmallScaleBias;
197 0 : itsStopAtLargeScaleNegative = other.itsStopAtLargeScaleNegative;
198 0 : itsStopPointMode = other.itsStopPointMode;
199 0 : itsDidStopPointMode = other.itsDidStopPointMode;
200 0 : itsJustStarting = other.itsJustStarting;
201 0 : itsStrengthOptimum = other.itsStrengthOptimum;
202 0 : itsTotalFlux=other.itsTotalFlux;
203 0 : itsMaskThreshold = other.itsMaskThreshold;
204 0 : psfShape_p.resize(0, false);
205 0 : psfShape_p=other.psfShape_p;
206 0 : noClean_p=other.noClean_p;
207 : }
208 0 : return *this;
209 : }
210 :
211 0 : MatrixCleaner::~MatrixCleaner()
212 : {
213 0 : destroyScales();
214 0 : if(!itsMask.null()) itsMask=0;
215 0 : }
216 :
217 :
218 0 : void MatrixCleaner::makeDirtyScales(){
219 :
220 0 : if(!itsScalesValid || itsNscales < 1 || itsDirty.null() || (itsNscales != Int(itsScaleXfrs.nelements())) )
221 0 : return;
222 :
223 0 : if( (psfShape_p) != (itsDirty->shape()))
224 0 : throw(AipsError("PSF and Dirty array are not of the same shape"));
225 : //No need of convolving for hogbom
226 0 : if(itsCleanType==CleanEnums::HOGBOM){
227 0 : itsDirtyConvScales.resize(1, true);
228 0 : itsDirtyConvScales[0]=Matrix<Float>(itsDirty->shape());
229 0 : itsDirtyConvScales[0]=*itsDirty;
230 0 : return;
231 : }
232 0 : Matrix<Complex> dirtyFT;
233 0 : FFTServer<Float,Complex> fft(itsDirty->shape());
234 0 : fft.fft0(dirtyFT, *itsDirty);
235 0 : itsDirtyConvScales.resize(itsNscales, true);
236 0 : Int scale=0;
237 : //Having problem with fftw and omp
238 : // #pragma omp parallel default(shared) private(scale) firstprivate(fft)
239 : {
240 : //#pragma omp for
241 : // Dirty*scale
242 0 : for (scale=0; scale<itsNscales; scale++) {
243 0 : Matrix<Complex> cWork;
244 : // Dirty * scale
245 : // cout << "scale " << scale << " itsScaleptr " << &(itsScaleXfrs[scale]) << "\n"<< endl;
246 :
247 0 : itsDirtyConvScales[scale]=Matrix<Float>(itsDirty->shape());
248 0 : cWork=((dirtyFT)*(itsScaleXfrs[scale]));
249 0 : fft.fft0((itsDirtyConvScales[scale]), cWork, false);
250 0 : fft.flip((itsDirtyConvScales[scale]), false, false);
251 0 : }
252 : }
253 :
254 0 : }
255 0 : void MatrixCleaner::update(const Matrix<Float> &dirty)
256 : {
257 0 : itsDirty->assign(dirty);
258 :
259 0 : LogIO os(LogOrigin("MatrixCleaner", "update()", WHERE));
260 :
261 :
262 :
263 : // Now we can redo the relevant convolutions
264 0 : makeDirtyScales();
265 0 : }
266 :
267 :
268 :
269 : // add a mask image
270 0 : void MatrixCleaner::setMask(Matrix<Float> & mask, const Float& maskThreshold)
271 : {
272 0 : itsMaskThreshold = maskThreshold;
273 0 : IPosition maskShape = mask.shape();
274 :
275 : //cerr << "Mask Shape " << mask.shape() << endl;
276 : // This is not needed after the first steps
277 0 : itsMask = new Matrix<Float>(mask.shape());
278 0 : itsMask->assign(mask);
279 0 : noClean_p=(max(*itsMask) < itsMaskThreshold) ? true : false;
280 :
281 :
282 :
283 0 : if (itsScalesValid) {
284 0 : makeScaleMasks();
285 : }
286 :
287 0 : }
288 :
289 0 : Bool MatrixCleaner::setcontrol(CleanEnums::CleanType cleanType,
290 : const Int niter,
291 : const Float gain,
292 : const Quantity& threshold)
293 : {
294 0 : return setcontrol(cleanType, niter, gain, threshold, Quantity(0.0, "%"));
295 : }
296 :
297 : // Set up the control parameters
298 0 : Bool MatrixCleaner::setcontrol(CleanEnums::CleanType cleanType,
299 : const Int niter,
300 : const Float gain,
301 : const Quantity& aThreshold,
302 : const Quantity& fThreshold)
303 : {
304 0 : itsCleanType=cleanType;
305 0 : itsMaxNiter=niter;
306 0 : itsGain=gain;
307 0 : itsThreshold=aThreshold;
308 0 : itsFracThreshold=fThreshold;
309 0 : return true;
310 : }
311 :
312 : // Set up speedup parameters
313 0 : void MatrixCleaner::speedup(const Float nDouble)
314 : {
315 0 : itsDoSpeedup=true;
316 0 : itsNDouble = nDouble;
317 0 : };
318 :
319 :
320 :
321 : // Do the clean as set up
322 0 : Int MatrixCleaner::clean(Matrix<Float>& model,
323 : Bool /*showProgress*/)
324 : {
325 0 : AlwaysAssert(model.shape()==itsDirty->shape(), AipsError);
326 :
327 0 : LogIO os(LogOrigin("MatrixCleaner", "clean()", WHERE));
328 :
329 0 : Float tmpMaximumResidual=0.0;
330 :
331 0 : Int nScalesToClean=itsNscales;
332 0 : if (itsCleanType==CleanEnums::HOGBOM) {
333 0 : os << LogIO::NORMAL1 << "Hogbom clean algorithm" << LogIO::POST;
334 0 : nScalesToClean=1;
335 : }
336 0 : else if (itsCleanType==CleanEnums::MULTISCALE) {
337 0 : if (nScalesToClean==1) {
338 0 : os << LogIO::NORMAL1 << "Multi-scale clean with only one scale" << LogIO::POST;
339 : }
340 : else {
341 0 : os << LogIO::NORMAL1 << "Multi-scale clean algorithm" << LogIO::POST;
342 : }
343 : }
344 :
345 : Int scale;
346 0 : Vector<Float> scaleBias(nScalesToClean);
347 0 : if (nScalesToClean > 1) {
348 0 : for (scale=0;scale<nScalesToClean;scale++) {
349 0 : scaleBias(scale) = 1 - itsSmallScaleBias *
350 0 : itsScaleSizes(scale)/itsScaleSizes(nScalesToClean-1);
351 0 : os << "scale " << scale+1 << " = " << itsScaleSizes(scale) << " pixels with bias = " << scaleBias(scale) << LogIO::POST;
352 : }
353 : } else {
354 0 : scaleBias(0) = 1.0;
355 : }
356 :
357 0 : AlwaysAssert(itsScalesValid, AipsError);
358 : ////no need to use all cores if possible
359 0 : Int nth=nScalesToClean;
360 : #ifdef _OPENMP
361 :
362 0 : nth=min(nth, omp_get_max_threads());
363 :
364 : #endif
365 : // Find the peaks of the convolved Psfs
366 0 : Vector<Float> maxPsfConvScales(nScalesToClean);
367 0 : Int naxes=model.shape().nelements();
368 0 : #pragma omp parallel default(shared) private(scale) num_threads(nth)
369 : {
370 : #pragma omp for
371 : for (scale=0;scale<nScalesToClean;scale++) {
372 : IPosition positionPeakPsfConvScales(naxes, 0);
373 :
374 : findMaxAbs(itsPsfConvScales[scale], maxPsfConvScales(scale),
375 : positionPeakPsfConvScales);
376 :
377 : // cout << "MAX " << scale << " " << positionPeakPsfConvScales
378 : // << " " << maxPsfConvScales(scale) << "\n"
379 : // << endl;
380 : }
381 : } //End pragma parallel
382 0 : for (scale=0;scale<nScalesToClean;scale++) {
383 0 : if ( maxPsfConvScales(scale) < 0.0) {
384 : os << "As Peak of PSF is negative, you should setscales again with a smaller scale size"
385 0 : << LogIO::SEVERE;
386 0 : return -1;
387 : }
388 : }
389 :
390 : // Define a subregion for the inner quarter
391 0 : IPosition blcDirty(model.shape().nelements(), 0);
392 0 : IPosition trcDirty(model.shape()-1);
393 :
394 0 : if(!itsMask.null()){
395 0 : os << "Cleaning using given mask" << LogIO::POST;
396 0 : if (itsMaskThreshold<0) {
397 : os << LogIO::NORMAL
398 : << "Mask thresholding is not used, values are interpreted as weights"
399 0 : <<LogIO::POST;
400 : } else {
401 : // a mask that does not allow for clean was sent
402 0 : if(noClean_p)
403 0 : return 0;
404 : os << LogIO::NORMAL
405 0 : << "Cleaning pixels with mask values above " << itsMaskThreshold
406 0 : << LogIO::POST;
407 : }
408 :
409 :
410 0 : Int nx=model.shape()(0);
411 0 : Int ny=model.shape()(1);
412 :
413 :
414 0 : AlwaysAssert(itsMask->shape()(0)==nx, AipsError);
415 0 : AlwaysAssert(itsMask->shape()(1)==ny, AipsError);
416 0 : Int xbeg=nx-1;
417 0 : Int ybeg=ny-1;
418 0 : Int xend=0;
419 0 : Int yend=0;
420 0 : for (Int iy=0;iy<ny;iy++) {
421 0 : for (Int ix=0;ix<nx;ix++) {
422 0 : if((*itsMask)(ix,iy)>0.000001) {
423 0 : xbeg=min(xbeg,ix);
424 0 : ybeg=min(ybeg,iy);
425 0 : xend=max(xend,ix);
426 0 : yend=max(yend,iy);
427 : }
428 : }
429 : }
430 :
431 0 : if (!itsIgnoreCenterBox) {
432 0 : if((xend - xbeg)>nx/2) {
433 0 : xbeg=nx/4-1; //if larger than quarter take inner of mask
434 0 : os << LogIO::WARN << "Mask span over more than half the x-axis: Considering inner half of the x-axis" << LogIO::POST;
435 : }
436 0 : if((yend - ybeg)>ny/2) {
437 0 : ybeg=ny/4-1;
438 0 : os << LogIO::WARN << "Mask span over more than half the y-axis: Considering inner half of the y-axis" << LogIO::POST;
439 : }
440 0 : xend=min(xend,xbeg+nx/2-1);
441 0 : yend=min(yend,ybeg+ny/2-1);
442 : }
443 :
444 : //blcDirty(0)=xbeg> 0 ? xbeg-1 : 0;
445 : //blcDirty(1)=ybeg > 0 ? ybeg-1 : 0;
446 0 : blcDirty(0)=xbeg;
447 0 : blcDirty(1)=ybeg;
448 0 : trcDirty(0)=xend;
449 0 : trcDirty(1)=yend;
450 : }
451 : else {
452 0 : if (itsIgnoreCenterBox) {
453 0 : os << LogIO::NORMAL << "Cleaning entire image" << LogIO::POST;
454 0 : os << LogIO::NORMAL1 << "as per MF/WF" << LogIO::POST; // ???
455 : }
456 : else {
457 0 : os << "Cleaning inner quarter of the image" << LogIO::POST;
458 0 : for (Int i=0;i<Int(model.shape().nelements());i++) {
459 0 : blcDirty(i)=model.shape()(i)/4;
460 0 : trcDirty(i)=blcDirty(i)+model.shape()(i)/2-1;
461 0 : if(trcDirty(i)<0) trcDirty(i)=1;
462 : }
463 : }
464 : }
465 0 : LCBox centerBox(blcDirty, trcDirty, model.shape());
466 :
467 0 : Block<Matrix<Float> > scaleMaskSubs;
468 0 : if (!itsMask.null()) {
469 0 : scaleMaskSubs.resize(itsNscales);
470 0 : for (Int is=0; is < itsNscales; is++) {
471 0 : scaleMaskSubs[is] = ((itsScaleMasks[is]))(blcDirty, trcDirty);
472 : }
473 : }
474 :
475 : // Start the iteration
476 0 : Vector<Float> maxima(nScalesToClean);
477 0 : Block<IPosition> posMaximum(nScalesToClean);
478 0 : Vector<Float> totalFluxScale(nScalesToClean); totalFluxScale=0.0;
479 0 : Float totalFlux=0.0;
480 0 : Int converged=0;
481 0 : Int stopPointModeCounter = 0;
482 0 : Int optimumScale=0;
483 0 : itsStrengthOptimum=0.0;
484 0 : IPosition positionOptimum(model.shape().nelements(), 0);
485 0 : os << "Starting iteration"<< LogIO::POST;
486 :
487 : //
488 0 : Int nx=model.shape()(0);
489 0 : Int ny=model.shape()(1);
490 0 : IPosition gip;
491 0 : gip = IPosition(2,nx,ny);
492 0 : casacore::Block<casacore::Matrix<casacore::Float> > vecWork_p;
493 0 : vecWork_p.resize(nScalesToClean);
494 :
495 0 : for(Int i=0;i<nScalesToClean;i++)
496 : {
497 0 : vecWork_p[i].resize(gip);
498 : }
499 : //
500 :
501 0 : itsIteration = itsStartingIter;
502 0 : for (Int ii=itsStartingIter; ii < itsMaxNiter; ii++) {
503 0 : itsIteration++;
504 : // Find the peak residual
505 0 : itsStrengthOptimum = 0.0;
506 0 : optimumScale = 0;
507 :
508 0 : #pragma omp parallel default(shared) private(scale) num_threads(nth)
509 : {
510 : #pragma omp for
511 : for (scale=0; scale<nScalesToClean; ++scale) {
512 : // Find absolute maximum for the dirty image
513 : // cout << "in omp loop for scale : " << scale << " : " << blcDirty << " : " << trcDirty << " :: " << itsDirtyConvScales.nelements() << endl;
514 : Matrix<Float> work = (vecWork_p[scale])(blcDirty,trcDirty);
515 : work = 0.0;
516 : work = work + (itsDirtyConvScales[scale])(blcDirty,trcDirty);
517 : maxima(scale)=0;
518 : posMaximum[scale]=IPosition(model.shape().nelements(), 0);
519 :
520 : if (!itsMask.null()) {
521 : findMaxAbsMask(vecWork_p[scale], itsScaleMasks[scale],
522 : maxima(scale), posMaximum[scale]);
523 : } else {
524 : findMaxAbs(vecWork_p[scale], maxima(scale), posMaximum[scale]);
525 : }
526 :
527 : // Remember to adjust the position for the window and for
528 : // the flux scale
529 : //cout << "scale " << scale << " maxPsfconvscale " << maxPsfConvScales(scale) << endl;
530 : //cout << "posmax " << posMaximum[scale] << " blcdir " << blcDirty << endl;
531 : maxima(scale)/=maxPsfConvScales(scale);
532 : maxima(scale) *= scaleBias(scale);
533 : maxima(scale) *= (itsDirtyConvScales[scale])(posMaximum[scale]); //makes maxima(scale) positive to ensure correct scale is selected in itsStrengthOptimum for loop (next for loop).
534 :
535 : //posMaximum[scale]+=blcDirty;
536 :
537 : }
538 : }//End parallel section
539 0 : for (scale=0; scale<nScalesToClean; scale++) {
540 0 : if(abs(maxima(scale))>abs(itsStrengthOptimum)) {
541 0 : optimumScale=scale;
542 0 : itsStrengthOptimum=maxima(scale);
543 0 : positionOptimum=posMaximum[scale];
544 : }
545 : }
546 :
547 : // These checks and messages are really only here as an early alert to us if SDAlgorithmBase::deconvolve does not skip all-0 images.
548 0 : if (abs(itsStrengthOptimum) == 0) {
549 0 : os << LogIO::SEVERE << "Attempting to divide 0. Scale maxima for scale " << optimumScale << "is 0." << LogIO::POST;
550 : } else {
551 0 : if (abs(scaleBias(optimumScale)) == 0) {
552 0 : os << LogIO::SEVERE << "Attempting to divide by 0. Scale bias for scale " << optimumScale << "is 0." << LogIO::POST;
553 : }
554 0 : itsStrengthOptimum /= scaleBias(optimumScale);
555 0 : itsStrengthOptimum /= (itsDirtyConvScales[optimumScale])(posMaximum[optimumScale]);
556 : }
557 :
558 0 : AlwaysAssert(optimumScale<nScalesToClean, AipsError);
559 :
560 : // Now add to the total flux
561 0 : totalFlux += (itsStrengthOptimum*itsGain);
562 0 : itsTotalFlux=totalFlux;
563 0 : totalFluxScale(optimumScale) += (itsStrengthOptimum*itsGain);
564 :
565 0 : if(ii==itsStartingIter ) {
566 0 : itsMaximumResidual=abs(itsStrengthOptimum);
567 0 : tmpMaximumResidual=itsMaximumResidual;
568 0 : os << "Initial maximum residual is " << itsMaximumResidual;
569 0 : if( !itsMask.null() ) { os << " within the mask "; }
570 0 : os << LogIO::POST;
571 : }
572 :
573 : // Various ways of stopping:
574 : // 1. stop if below threshold
575 0 : if(abs(itsStrengthOptimum)<threshold() ) {
576 0 : os << "Reached stopping threshold " << threshold() << " at iteration "<<
577 0 : ii << LogIO::POST;
578 0 : os << "Optimum flux is " << abs(itsStrengthOptimum) << LogIO::POST;
579 0 : converged = 1;
580 0 : break;
581 : }
582 : // 2. negatives on largest scale?
583 0 : if ((nScalesToClean > 1) && itsStopAtLargeScaleNegative &&
584 0 : optimumScale == (nScalesToClean-1) &&
585 0 : itsStrengthOptimum < 0.0) {
586 0 : os << "Reached negative on largest scale" << LogIO::POST;
587 0 : converged = -2;
588 0 : break;
589 : }
590 : // 3. stop point mode at work
591 0 : if (itsStopPointMode > 0) {
592 0 : if (optimumScale == 0) {
593 0 : stopPointModeCounter++;
594 : } else {
595 0 : stopPointModeCounter = 0;
596 : }
597 0 : if (stopPointModeCounter >= itsStopPointMode) {
598 : os << "Cleaned " << stopPointModeCounter <<
599 : " consecutive components from the smallest scale, stopping prematurely"
600 0 : << LogIO::POST;
601 0 : itsDidStopPointMode = true;
602 0 : converged = -1;
603 0 : break;
604 : }
605 : }
606 : //4. Diverging large scale
607 : //If actual value is 50% above the maximum residual. ..good chance it will not recover at this stage
608 0 : if(((abs(itsStrengthOptimum)-abs(tmpMaximumResidual)) > (abs(tmpMaximumResidual)/2.0))
609 0 : && !(itsStopAtLargeScaleNegative)){
610 : os << "Diverging due to large scale?"
611 0 : << LogIO::POST;
612 : //clean is diverging most probably due to the large scale
613 0 : converged=-2;
614 0 : break;
615 : }
616 : //5. Diverging for some other reason; may just need another CS-style reconciling
617 0 : if((abs(itsStrengthOptimum)-abs(tmpMaximumResidual)) > (abs(tmpMaximumResidual)/2.0)){
618 : os << "Diverging due to unknown reason"
619 0 : << LogIO::POST;
620 0 : converged=-3;
621 0 : break;
622 : }
623 :
624 :
625 : /*
626 : if(progress) {
627 : progress->info(false, itsIteration, itsMaxNiter, maxima,
628 : posMaximum, itsStrengthOptimum,
629 : optimumScale, positionOptimum,
630 : totalFlux, totalFluxScale,
631 : itsJustStarting );
632 : itsJustStarting = false;
633 : } else*/ {
634 0 : if (itsIteration == itsStartingIter + 1) {
635 0 : os << "iteration MaximumResidual CleanedFlux" << LogIO::POST;
636 : }
637 0 : if ((itsIteration % (itsMaxNiter/10 > 0 ? itsMaxNiter/10 : 1)) == 0) {
638 : //Good place to re-up the fiducial maximum residual
639 : //tmpMaximumResidual=abs(itsStrengthOptimum);
640 0 : os << itsIteration <<" "<<itsStrengthOptimum<<" "
641 0 : << totalFlux <<LogIO::POST;
642 : }
643 : }
644 :
645 : Float scaleFactor;
646 0 : scaleFactor=itsGain*itsStrengthOptimum;
647 :
648 : // Continuing: subtract the peak that we found from all dirty images
649 : // Define a subregion so that that the peak is centered
650 0 : IPosition support(model.shape());
651 0 : support(0)=max(Int(itsScaleSizes(itsNscales-1)+0.5), support(0));
652 0 : support(1)=max(Int(itsScaleSizes(itsNscales-1)+0.5), support(1));
653 :
654 0 : IPosition inc(model.shape().nelements(), 1);
655 : //cout << "support " << support.asVector() << endl;
656 : //support(0)=1024;
657 : //support(1)=1024;
658 : //support(0)=min(Int(support(0)), Int(trcDirty(0)-blcDirty(0)));
659 : //support(1)=min(Int(support(1)), Int(trcDirty(1)-blcDirty(1)));
660 : // support(0)=min(Int(support(0)), (trcDirty(0)-blcDirty(0)+
661 : // Int(2*abs(positionOptimum(0)-blcDirty(0)/2.0-trcDirty(0)/2.0))));
662 : //support(1)=min(Int(support(1)), (trcDirty(1)-blcDirty(1)+
663 : // Int(2*abs(positionOptimum(1)-blcDirty(1)/2.0-trcDirty(1)/2.0))));
664 :
665 0 : IPosition blc(positionOptimum-support/2);
666 0 : IPosition trc(positionOptimum+support/2-1);
667 0 : LCBox::verify(blc, trc, inc, model.shape());
668 :
669 : //cout << "blc " << blc.asVector() << " trc " << trc.asVector() << endl;
670 :
671 0 : IPosition blcPsf(blc+itsPositionPeakPsf-positionOptimum);
672 0 : IPosition trcPsf(trc+itsPositionPeakPsf-positionOptimum);
673 0 : LCBox::verify(blcPsf, trcPsf, inc, model.shape());
674 0 : makeBoxesSameSize(blc,trc,blcPsf,trcPsf);
675 : // cout << "blcPsf " << blcPsf.asVector() << " trcPsf " << trcPsf.asVector() << endl;
676 : //cout << "blc " << blc.asVector() << " trc " << trc.asVector() << endl;
677 : // LCBox subRegion(blc, trc, model.shape());
678 : // LCBox subRegionPsf(blcPsf, trcPsf, model.shape());
679 :
680 0 : Matrix<Float> modelSub=model(blc, trc);
681 0 : Matrix<Float> scaleSub=(itsScales[optimumScale])(blcPsf,trcPsf);
682 :
683 :
684 : // Now do the addition of this scale to the model image....
685 0 : modelSub += scaleFactor*scaleSub;
686 :
687 0 : #pragma omp parallel default(shared) private(scale) num_threads(nth)
688 : {
689 : #pragma omp for
690 : for (scale=0;scale<nScalesToClean; ++scale) {
691 :
692 : Matrix<Float> dirtySub=(itsDirtyConvScales[scale])(blc,trc);
693 : //AlwaysAssert(itsPsfConvScales[index(scale,optimumScale)], AipsError);
694 : Matrix<Float> psfSub=(itsPsfConvScales[index(scale,optimumScale)])(blcPsf, trcPsf);
695 : dirtySub -= scaleFactor*psfSub;
696 :
697 : }
698 : }//End parallel
699 :
700 0 : blcDirty = blc;
701 0 : trcDirty = trc;
702 0 : }
703 : // End of iteration
704 :
705 0 : for (scale=0;scale<nScalesToClean;scale++) {
706 : os << LogIO::NORMAL
707 0 : << " " << scale << " " << totalFluxScale(scale)
708 0 : << LogIO::POST;
709 : }
710 : // Finish off the plot, etc.
711 : /*
712 : if(progress) {
713 : progress->info(true, itsIteration, itsMaxNiter, maxima, posMaximum,
714 : itsStrengthOptimum,
715 : optimumScale, positionOptimum,
716 : totalFlux, totalFluxScale);
717 : }
718 : */
719 :
720 0 : if(!converged) {
721 0 : os << "Failed to reach stopping threshold" << LogIO::POST;
722 : }
723 :
724 0 : return converged;
725 0 : }
726 :
727 0 : Bool MatrixCleaner::findPSFMaxAbs(const Matrix<Float>& lattice,
728 : Float& maxAbs,
729 : IPosition& posMaxAbs,
730 : const Int& supportSize)
731 : {
732 0 : LogIO os(LogOrigin("MatrixCleaner", "findPSFMaxAbs()", WHERE));
733 0 : posMaxAbs = IPosition(lattice.shape().nelements(), 0);
734 0 : maxAbs=0.0;
735 :
736 0 : Float maxVal=0.0;
737 0 : IPosition psf2DShape(lattice.shape());
738 0 : Int blc0 = (psf2DShape(0) > supportSize) ? psf2DShape(0)/2 - supportSize/2 : 0;
739 0 : Int blc1 = (psf2DShape(1) > supportSize) ? psf2DShape(1)/2 - supportSize/2 : 0;
740 0 : Int trc0 = (psf2DShape(0) > supportSize) ? psf2DShape(0)/2 + supportSize/2 : psf2DShape(0)-1;
741 :
742 0 : Int trc1 = (psf2DShape(1) > supportSize) ? (psf2DShape(1)/2 + supportSize/2) : (psf2DShape(1)-1) ;
743 :
744 :
745 : // cerr << "####### " << blc0 << " " << blc1 << " " << trc0 << " " << trc1 << endl;
746 : // cerr << "Max of lattice " << max(lattice) << " min " << min(lattice) << endl;
747 0 : for (Int j=blc1; j < trc1; ++j)
748 0 : for (Int i=blc0 ; i < trc0; ++i)
749 0 : if ((maxAbs = abs(lattice(i,j))) > maxVal)
750 : {
751 0 : maxVal = maxAbs;
752 0 : posMaxAbs(0)=i; posMaxAbs(1)=j;
753 : }
754 0 : maxAbs=maxVal;
755 : //cerr << "######## " << posMaxAbs << " " << maxVal << endl;
756 0 : return true;
757 0 : }
758 :
759 0 : Bool MatrixCleaner::findMaxAbs(const Matrix<Float>& lattice,
760 : Float& maxAbs,
761 : IPosition& posMaxAbs)
762 : {
763 :
764 0 : posMaxAbs = IPosition(lattice.shape().nelements(), 0);
765 0 : maxAbs=0.0;
766 :
767 : Float minVal;
768 0 : IPosition posmin(lattice.shape().nelements(), 0);
769 0 : minMax(minVal, maxAbs, posmin, posMaxAbs, lattice);
770 : //cout << "min " << minVal << " " << maxAbs << " " << max(lattice) << endl;
771 0 : if(abs(minVal) > abs(maxAbs)){
772 0 : maxAbs=minVal;
773 0 : posMaxAbs=posmin;
774 : }
775 0 : return true;
776 0 : }
777 :
778 :
779 :
780 :
781 :
782 0 : Bool MatrixCleaner::findMaxAbsMask(const Matrix<Float>& lattice,
783 : const Matrix<Float>& mask,
784 : Float& maxAbs,
785 : IPosition& posMaxAbs)
786 : {
787 :
788 0 : posMaxAbs = IPosition(lattice.shape().nelements(), 0);
789 0 : maxAbs=0.0;
790 : Float minVal;
791 0 : IPosition posmin(lattice.shape().nelements(), 0);
792 0 : minMaxMasked(minVal, maxAbs, posmin, posMaxAbs, lattice, mask);
793 0 : if(abs(minVal) > abs(maxAbs)){
794 0 : maxAbs=minVal;
795 0 : posMaxAbs=posmin;
796 : }
797 :
798 0 : return true;
799 0 : }
800 :
801 :
802 :
803 0 : Bool MatrixCleaner::setscales(const Int nscales, const Float scaleInc)
804 : {
805 0 : LogIO os(LogOrigin("deconvolver", "setscales()", WHERE));
806 :
807 :
808 0 : itsNscales=nscales;
809 0 : if(itsNscales<1) {
810 0 : os << "Using default of 5 scales" << LogIO::POST;
811 0 : itsNscales=5;
812 : }
813 :
814 0 : Vector<Float> scaleSizes(itsNscales);
815 :
816 : // Validate scales
817 0 : os << "Creating " << itsNscales << " scales" << LogIO::POST;
818 0 : scaleSizes(0) = 0.00001 * scaleInc;
819 0 : os << "scale 1 = 0.0 arcsec" << LogIO::POST;
820 0 : for (Int scale=1; scale<itsNscales;scale++) {
821 0 : scaleSizes(scale) =
822 0 : scaleInc * pow(10.0, (Float(scale)-2.0)/2.0);
823 0 : os << "scale " << scale+1 << " = " << scaleSizes(scale)
824 0 : << " arcsec" << LogIO::POST;
825 : }
826 :
827 0 : return setscales(scaleSizes);
828 :
829 0 : }
830 :
831 :
832 0 : void MatrixCleaner::setDirty(const Matrix<Float>& dirty){
833 :
834 0 : itsDirty=new Matrix<Float>(dirty.shape());
835 0 : itsDirty->assign(dirty);
836 :
837 :
838 :
839 0 : }
840 :
841 : //Define the scales without doing anything else
842 : // user will call make makePsfScales and makeDirtyScales like an adult in the know
843 :
844 0 : void MatrixCleaner::defineScales(const Vector<Float>& scaleSizes){
845 0 : if(itsScales.nelements()>0) {
846 0 : destroyScales();
847 : }
848 :
849 0 : destroyMasks();
850 0 : itsNscales=scaleSizes.nelements();
851 0 : itsScaleSizes.resize(itsNscales);
852 0 : itsScaleSizes=scaleSizes; // make a copy that we can call our own
853 0 : GenSort<Float>::sort(itsScaleSizes);
854 0 : itsScalesValid=false;
855 0 : }
856 :
857 0 : void MatrixCleaner::makePsfScales(){
858 0 : LogIO os(LogOrigin("MatrixCleaner", "makePsfScales()", WHERE));
859 0 : if(itsNscales < 1)
860 0 : throw(AipsError("Scales have to be set"));
861 0 : if(itsXfr.null())
862 0 : throw(AipsError("Psf is not defined"));
863 0 : destroyScales();
864 0 : itsScales.resize(itsNscales, true);
865 0 : itsScaleXfrs.resize(itsNscales, true);
866 0 : itsPsfConvScales.resize((itsNscales+1)*(itsNscales+1), true);
867 0 : FFTServer<Float,Complex> fft(psfShape_p);
868 0 : Int scale=0;
869 0 : for(scale=0; scale<itsNscales;scale++) {
870 0 : itsScales[scale] = Matrix<Float>(psfShape_p);
871 0 : makeScale(itsScales[scale], itsScaleSizes(scale));
872 0 : itsScaleXfrs[scale] = Matrix<Complex> ();
873 0 : fft.fft0(itsScaleXfrs[scale], itsScales[scale]);
874 : }
875 0 : Matrix<Complex> cWork;
876 :
877 0 : for (scale=0; scale<itsNscales;scale++) {
878 0 : os << "Calculating convolutions for scale " << scale << LogIO::POST;
879 : //PSF * scale
880 0 : itsPsfConvScales[scale] = Matrix<Float>(psfShape_p);
881 0 : cWork=((*itsXfr)*(itsScaleXfrs[scale])*(itsScaleXfrs[scale]));
882 : //cout << "shape " << cWork.shape() << " " << itsPsfConvScales[scale].shape() << endl;
883 :
884 0 : fft.fft0((itsPsfConvScales[scale]), cWork, false);
885 0 : fft.flip(itsPsfConvScales[scale], false, false);
886 :
887 : //cout << "psf scale " << scale << " " << max(itsPsfConvScales[scale]) << " " << min(itsPsfConvScales[scale]) << endl;
888 :
889 0 : for (Int otherscale=scale;otherscale<itsNscales;otherscale++) {
890 :
891 0 : AlwaysAssert(index(scale, otherscale)<Int(itsPsfConvScales.nelements()),
892 : AipsError);
893 :
894 : // PSF * scale * otherscale
895 0 : itsPsfConvScales[index(scale,otherscale)] =Matrix<Float>(psfShape_p);
896 0 : cWork=((*itsXfr)*(itsScaleXfrs[scale])*(itsScaleXfrs[otherscale]));
897 0 : fft.fft0(itsPsfConvScales[index(scale,otherscale)], cWork, false);
898 : //For some reason this complex->real fft does not need a flip ...may be because conj(a)*a is real
899 : //fft.flip(*itsPsfConvScales[index(scale,otherscale)], false, false);
900 : }
901 : }
902 :
903 0 : itsScalesValid=true;
904 :
905 0 : }
906 :
907 : // We calculate all the scales and the corresponding convolutions
908 : // and cross convolutions.
909 :
910 0 : Bool MatrixCleaner::setscales(const Vector<Float>& scaleSizes)
911 : {
912 0 : LogIO os(LogOrigin("MatrixCleaner", "setscales()", WHERE));
913 :
914 : Int scale;
915 :
916 0 : defineScales(scaleSizes);
917 :
918 : // Residual, psf, and mask, plus cross terms
919 : // e.g. for 5 scales this is 45. for 6 it is 60.
920 0 : Int nImages=3*itsNscales+itsNscales*(itsNscales+1);
921 0 : os << "Expect to use " << nImages << " scratch images" << LogIO::POST;
922 :
923 : // Now we can update the size of memory allocated
924 0 : itsMemoryMB=0.5*Double(HostInfo::memoryTotal()/1024)/Double(nImages);
925 0 : os << "Maximum memory allocated per image " << itsMemoryMB << "MB" << LogIO::POST;
926 :
927 0 : itsDirtyConvScales.resize(itsNscales);
928 0 : itsScaleMasks.resize(itsNscales);
929 0 : itsScaleXfrs.resize(itsNscales);
930 0 : itsPsfConvScales.resize((itsNscales+1)*(itsNscales+1));
931 0 : for(scale=0; scale<itsNscales;scale++) {
932 0 : itsDirtyConvScales[scale].resize();
933 0 : itsScaleMasks[scale].resize();
934 0 : itsScaleXfrs[scale].resize();
935 : }
936 0 : for(scale=0; scale<((itsNscales+1)*(itsNscales+1));scale++) {
937 0 : itsPsfConvScales[scale].resize();
938 : }
939 :
940 0 : AlwaysAssert(!itsDirty.null(), AipsError);
941 :
942 0 : FFTServer<Float,Complex> fft(itsDirty->shape());
943 :
944 0 : Matrix<Complex> dirtyFT;
945 0 : fft.fft0(dirtyFT, *itsDirty);
946 :
947 : ///Having problem with fftw with openmp
948 : //#pragma parallel default(shared) private(scale) firstprivate(fft)
949 : {
950 : //#pragma omp for
951 0 : for (scale=0; scale<itsNscales;scale++) {
952 0 : os << "Calculating scale image and Fourier transform for scale " << scale << LogIO::POST;
953 : //cout << "Calculating scale image and Fourier transform for scale " << scale << endl;
954 0 : itsScales[scale] = Matrix<Float>(itsDirty->shape());
955 : //AlwaysAssert(itsScales[scale], AipsError);
956 : // First make the scale
957 0 : makeScale(itsScales[scale], scaleSizes(scale));
958 0 : itsScaleXfrs[scale] = Matrix<Complex> ();
959 0 : fft.fft0(itsScaleXfrs[scale], itsScales[scale]);
960 : }
961 : }
962 :
963 : // Now we can do all the convolutions
964 0 : Matrix<Complex> cWork;
965 0 : for (scale=0; scale<itsNscales;scale++) {
966 0 : os << "Calculating convolutions for scale " << scale << LogIO::POST;
967 :
968 : // PSF * scale
969 0 : itsPsfConvScales[scale] = Matrix<Float>(itsDirty->shape());
970 :
971 0 : cWork=((*itsXfr)*(itsScaleXfrs[scale]));
972 :
973 0 : fft.fft0((itsPsfConvScales[scale]), cWork, false);
974 0 : fft.flip(itsPsfConvScales[scale], false, false);
975 :
976 0 : itsDirtyConvScales[scale] = Matrix<Float>(itsDirty->shape());
977 0 : cWork=((dirtyFT)*(itsScaleXfrs[scale]));
978 0 : fft.fft0(itsDirtyConvScales[scale], cWork, false);
979 0 : fft.flip(itsDirtyConvScales[scale], false, false);
980 : ///////////
981 : /*
982 : {
983 : String axisName = "TabularDoggies1";
984 : String axisUnit = "km";
985 : Double crval = 10.12;
986 : Double crpix = -128.32;
987 : Double cdelt = 3.145;
988 : TabularCoordinate tab1(crval, cdelt, crpix, axisUnit, axisName);
989 : TabularCoordinate tab2(crval, cdelt, crpix, axisUnit, "Dogma");
990 : CoordinateSystem csys;
991 : csys.addCoordinate(tab1);
992 : csys.addCoordinate(tab2);
993 : PagedImage<Float> limage(itsPsfConvScales[scale].shape(), csys, "psfconvscale_"+String::toString(scale));
994 : limage.put(itsPsfConvScales[scale]);
995 : }
996 : */
997 : ///////////
998 :
999 0 : for (Int otherscale=scale;otherscale<itsNscales;otherscale++) {
1000 :
1001 0 : AlwaysAssert(index(scale, otherscale)<Int(itsPsfConvScales.nelements()),
1002 : AipsError);
1003 :
1004 : // PSF * scale * otherscale
1005 0 : itsPsfConvScales[index(scale,otherscale)] =Matrix<Float>(itsDirty->shape());
1006 0 : cWork=((*itsXfr)*conj(itsScaleXfrs[scale])*(itsScaleXfrs[otherscale]));
1007 0 : fft.fft0(itsPsfConvScales[index(scale,otherscale)], cWork, false);
1008 : //fft.flip(*itsPsfConvScales[index(scale,otherscale)], false, false);
1009 : }
1010 : }
1011 :
1012 0 : itsScalesValid=true;
1013 :
1014 0 : if (!itsMask.null()) {
1015 0 : makeScaleMasks();
1016 : }
1017 :
1018 0 : return true;
1019 0 : }
1020 :
1021 : // Make a single scale size image
1022 0 : void MatrixCleaner::makeScale(Matrix<Float>& iscale, const Float& scaleSize)
1023 : {
1024 :
1025 0 : Int nx=iscale.shape()(0);
1026 0 : Int ny=iscale.shape()(1);
1027 : //Matrix<Float> iscale(nx, ny);
1028 0 : iscale=0.0;
1029 :
1030 0 : Double refi=nx/2;
1031 0 : Double refj=ny/2;
1032 :
1033 0 : if(scaleSize==0.0) {
1034 0 : iscale(Int(refi), Int(refj)) = 1.0;
1035 : }
1036 : else {
1037 0 : AlwaysAssert(scaleSize>0.0,AipsError);
1038 :
1039 0 : Int mini = max( 0, (Int)(refi-scaleSize));
1040 0 : Int maxi = min(nx-1, (Int)(refi+scaleSize));
1041 0 : Int minj = max( 0, (Int)(refj-scaleSize));
1042 0 : Int maxj = min(ny-1, (Int)(refj+scaleSize));
1043 :
1044 0 : Float ypart=0.0;
1045 0 : Float volume=0.0;
1046 0 : Float rad2=0.0;
1047 0 : Float rad=0.0;
1048 :
1049 0 : for (Int j=minj;j<=maxj;j++) {
1050 0 : ypart = square( (refj - (Double)(j)) / scaleSize );
1051 0 : for (Int i=mini;i<=maxi;i++) {
1052 0 : rad2 = ypart + square( (refi - (Double)(i)) / scaleSize );
1053 0 : if (rad2 < 1.0) {
1054 0 : if (rad2 <= 0.0) {
1055 0 : rad = 0.0;
1056 : } else {
1057 0 : rad = sqrt(rad2);
1058 : }
1059 0 : iscale(i,j) = (1.0 - rad2) * spheroidal(rad);
1060 0 : volume += iscale(i,j);
1061 : } else {
1062 0 : iscale(i,j) = 0.0;
1063 : }
1064 : }
1065 : }
1066 0 : iscale/=volume;
1067 : }
1068 : //scale.putSlice(iscale, IPosition(scale.ndim(),0), IPosition(scale.ndim(),1));
1069 0 : }
1070 :
1071 : // Calculate the spheroidal function
1072 0 : Float MatrixCleaner::spheroidal(Float nu) {
1073 :
1074 0 : if (nu <= 0) {
1075 0 : return 1.0;
1076 0 : } else if (nu >= 1.0) {
1077 0 : return 0.0;
1078 : } else {
1079 0 : uInt np = 5;
1080 0 : uInt nq = 3;
1081 0 : Matrix<Float> p(np, 2);
1082 0 : Matrix<Float> q(nq, 2);
1083 0 : p(0,0) = 8.203343e-2;
1084 0 : p(1,0) = -3.644705e-1;
1085 0 : p(2,0) = 6.278660e-1;
1086 0 : p(3,0) = -5.335581e-1;
1087 0 : p(4,0) = 2.312756e-1;
1088 0 : p(0,1) = 4.028559e-3;
1089 0 : p(1,1) = -3.697768e-2;
1090 0 : p(2,1) = 1.021332e-1;
1091 0 : p(3,1) = -1.201436e-1;
1092 0 : p(4,1) = 6.412774e-2;
1093 0 : q(0,0) = 1.0000000e0;
1094 0 : q(1,0) = 8.212018e-1;
1095 0 : q(2,0) = 2.078043e-1;
1096 0 : q(0,1) = 1.0000000e0;
1097 0 : q(1,1) = 9.599102e-1;
1098 0 : q(2,1) = 2.918724e-1;
1099 0 : uInt part = 0;
1100 0 : Float nuend = 0.0;
1101 0 : if (nu >= 0.0 && nu < 0.75) {
1102 0 : part = 0;
1103 0 : nuend = 0.75;
1104 0 : } else if (nu >= 0.75 && nu <= 1.00) {
1105 0 : part = 1;
1106 0 : nuend = 1.0;
1107 : }
1108 :
1109 0 : Float top = p(0,part);
1110 0 : Float delnusq = pow(nu,2.0) - pow(nuend,2.0);
1111 : uInt k;
1112 0 : for (k=1; k<np; k++) {
1113 0 : top += p(k, part) * pow(delnusq, (Float)k);
1114 : }
1115 0 : Float bot = q(0, part);
1116 0 : for (k=1; k<nq; k++) {
1117 0 : bot += q(k,part) * pow(delnusq, (Float)k);
1118 : }
1119 :
1120 0 : if (bot != 0.0) {
1121 0 : return (top/bot);
1122 : } else {
1123 0 : return 0.0;
1124 : }
1125 0 : }
1126 : }
1127 :
1128 :
1129 : // Calculate index into PsfConvScales
1130 0 : Int MatrixCleaner::index(const Int scale, const Int otherscale) {
1131 0 : if(otherscale>scale) {
1132 0 : return scale + itsNscales*(otherscale+1);
1133 : }
1134 : else {
1135 0 : return otherscale + itsNscales*(scale+1);
1136 : }
1137 : }
1138 :
1139 :
1140 0 : Bool MatrixCleaner::destroyScales()
1141 : {
1142 0 : if(!itsScalesValid) return true;
1143 0 : for(uInt scale=0; scale<itsScales.nelements();scale++) {
1144 : //if(itsScales[scale]) delete itsScales[scale];
1145 0 : itsScales[scale].resize();
1146 : }
1147 0 : for(uInt scale=0; scale<itsScaleXfrs.nelements();scale++) {
1148 : //if(itsScaleXfrs[scale]) delete itsScaleXfrs[scale];
1149 0 : itsScaleXfrs[scale].resize();
1150 : }
1151 0 : for(uInt scale=0; scale<itsDirtyConvScales.nelements();scale++) {
1152 : //if(itsDirtyConvScales[scale]) delete itsDirtyConvScales[scale];
1153 0 : itsDirtyConvScales[scale].resize();
1154 : }
1155 0 : for(uInt scale=0; scale<itsPsfConvScales.nelements();scale++) {
1156 : //if(itsPsfConvScales[scale]) delete itsPsfConvScales[scale];
1157 0 : itsPsfConvScales[scale].resize();
1158 : }
1159 0 : destroyMasks();
1160 0 : itsScales.resize(0, true);
1161 0 : itsDirtyConvScales.resize(0,true);
1162 0 : itsPsfConvScales.resize(0, true);
1163 0 : itsScalesValid=false;
1164 0 : return true;
1165 : }
1166 :
1167 :
1168 :
1169 0 : Bool MatrixCleaner::destroyMasks()
1170 : {
1171 0 : for(uInt scale=0; scale<itsScaleMasks.nelements();scale++) {
1172 : //if(itsScaleMasks[scale]) delete itsScaleMasks[scale];
1173 0 : itsScaleMasks[scale].resize();
1174 : }
1175 0 : itsScaleMasks.resize(0);
1176 0 : return true;
1177 : };
1178 :
1179 0 : void MatrixCleaner::unsetMask()
1180 : {
1181 0 : destroyMasks();
1182 0 : if(!itsMask.null())
1183 0 : itsMask=0;
1184 0 : noClean_p=false;
1185 0 : }
1186 :
1187 :
1188 :
1189 : // Set up the masks for the various scales
1190 : // This really only works for well behaved (ie, non-concave) masks.
1191 : // with only 1.0 or 0.0 values, and assuming the Scale images have
1192 : // a finite extent equal to +/- itsScaleSizes(scale)
1193 :
1194 0 : Bool MatrixCleaner::makeScaleMasks()
1195 : {
1196 0 : LogIO os(LogOrigin("MatrixCleaner", "makeScaleMasks()", WHERE));
1197 : Int scale;
1198 :
1199 :
1200 0 : if(!itsScalesValid) {
1201 : os << "Scales are not yet set - cannot set scale masks"
1202 0 : << LogIO::EXCEPTION;
1203 : }
1204 :
1205 :
1206 0 : destroyMasks();
1207 :
1208 0 : if(itsMask.null() || noClean_p)
1209 0 : return false;
1210 :
1211 0 : AlwaysAssert((itsMask->shape() == psfShape_p), AipsError);
1212 : //cout << "itsmask " << itsMask->shape() << endl;
1213 0 : FFTServer<Float,Complex> fft(itsMask->shape());
1214 :
1215 0 : Matrix<Complex> maskFT;
1216 0 : fft.fft0(maskFT, *itsMask);
1217 0 : itsScaleMasks.resize(itsNscales, true);
1218 : // Now we can do all the convolutions
1219 0 : Matrix<Complex> cWork;
1220 0 : for (scale=0; scale<itsNscales;scale++) {
1221 : //AlwaysAssert(itsScaleXfrs[scale], AipsError);
1222 0 : os << "Calculating mask convolution for scale " << scale << LogIO::POST;
1223 :
1224 : // Mask * scale
1225 : // Allow only 10% overlap by default, hence 0.9 is a default mask threshold
1226 : // if thresholding is not used, just extract the real part of the complex mask
1227 0 : itsScaleMasks[scale] = Matrix<Float>(itsMask->shape());
1228 0 : cWork=((maskFT)*(itsScaleXfrs[scale]));
1229 0 : fft.fft0(itsScaleMasks[scale], cWork, false);
1230 0 : fft.flip(itsScaleMasks[scale], false, false);
1231 0 : for (Int j=0 ; j < (itsMask->shape())(1); ++j){
1232 0 : for (Int k =0 ; k < (itsMask->shape())(0); ++k){
1233 0 : if(itsMaskThreshold > 0){
1234 0 : (itsScaleMasks[scale])(k,j) = (itsScaleMasks[scale])(k,j) > itsMaskThreshold ? 1.0 : 0.0;
1235 : }
1236 : }
1237 : }
1238 0 : Float mysum = sum(itsScaleMasks[scale] );
1239 0 : if (mysum <= 0.1) {
1240 : os << LogIO::WARN << "Ignoring scale " << scale <<
1241 0 : " since it is too large to fit within the mask" << LogIO::POST;
1242 : }
1243 :
1244 : }
1245 :
1246 0 : Int nx=itsScaleMasks[0].shape()(0);
1247 0 : Int ny=itsScaleMasks[0].shape()(1);
1248 :
1249 : /* Set the edges of the masks according to the scale size */
1250 : // Set the values OUTSIDE the box to zero....
1251 0 : for(Int scale=0;scale<itsNscales;scale++)
1252 : {
1253 0 : Int border = (Int)(itsScaleSizes[scale]*1.5);
1254 : // bottom
1255 0 : IPosition blc1(2, 0 , 0 );
1256 0 : IPosition trc1(2,nx-1, border );
1257 0 : IPosition inc1(2, 1);
1258 0 : LCBox::verify(blc1,trc1,inc1,itsScaleMasks[scale].shape());
1259 0 : (itsScaleMasks[scale])(blc1,trc1) = 0.0;
1260 : // top
1261 0 : blc1[0]=0; blc1[1]=ny-border-1;
1262 0 : trc1[0]=nx-1; trc1[1]=ny-1;
1263 0 : LCBox::verify(blc1,trc1,inc1,itsScaleMasks[scale].shape());
1264 0 : (itsScaleMasks[scale])(blc1,trc1) = 0.0;
1265 : // left
1266 0 : blc1[0]=0; blc1[1]=border;
1267 0 : trc1[0]=border; trc1[1]=ny-border-1;
1268 0 : LCBox::verify(blc1,trc1,inc1,itsScaleMasks[scale].shape());
1269 0 : (itsScaleMasks[scale])(blc1,trc1) = 0.0;
1270 : // right
1271 0 : blc1[0]=nx-border-1; blc1[1]=border;
1272 0 : trc1[0]=nx; trc1[1]=ny-border-1;
1273 0 : LCBox::verify(blc1,trc1,inc1,itsScaleMasks[scale].shape());
1274 0 : (itsScaleMasks[scale])(blc1,trc1) = 0.0;
1275 0 : }
1276 :
1277 0 : return true;
1278 0 : }
1279 :
1280 :
1281 0 : Matrix<casacore::Float> MatrixCleaner::residual(const Matrix<Float>& model){
1282 0 : FFTServer<Float,Complex> fft(model.shape());
1283 0 : if(!itsDirty.null() && (model.shape() != (itsDirty->shape())))
1284 0 : throw(AipsError("MatrixCleaner: model size has to be consistent with dirty image size"));
1285 0 : if(itsXfr.null())
1286 0 : throw(AipsError("MatrixCleaner: need to know the psf to calculate the residual"));
1287 :
1288 0 : Matrix<Complex> modelFT;
1289 : ///Convolve model with psf
1290 0 : fft.fft0(modelFT, model);
1291 0 : modelFT= (*itsXfr) * modelFT;
1292 0 : Matrix<Float> newResidual(itsDirty->shape());
1293 0 : fft.fft0(newResidual, modelFT, false);
1294 0 : fft.flip(newResidual, false, false);
1295 0 : newResidual=(*itsDirty)-newResidual;
1296 0 : return newResidual;
1297 0 : }
1298 :
1299 :
1300 0 : Float MatrixCleaner::threshold() const
1301 : {
1302 0 : if (! itsDoSpeedup) {
1303 0 : return max(itsFracThreshold.get("%").getValue() * itsMaximumResidual /100.0,
1304 0 : itsThreshold.get("Jy").getValue());
1305 : } else {
1306 0 : const Float factor = exp( (Float)( itsIteration - itsStartingIter )/ itsNDouble )
1307 0 : / 2.7182818;
1308 0 : return factor * max(itsFracThreshold.get("%").getValue() * itsMaximumResidual /100.0,
1309 0 : itsThreshold.get("Jy").getValue());
1310 : }
1311 : }
1312 :
1313 :
1314 :
1315 0 : void MatrixCleaner::makeBoxesSameSize(IPosition& blc1, IPosition& trc1,
1316 : IPosition &blc2, IPosition& trc2)
1317 : {
1318 0 : const IPosition shape1 = trc1 - blc1;
1319 0 : const IPosition shape2 = trc2 - blc2;
1320 :
1321 0 : AlwaysAssert(shape1.nelements() == shape2.nelements(), AipsError);
1322 :
1323 0 : if (shape1 == shape2) {
1324 0 : return;
1325 : }
1326 0 : for (uInt i=0;i<shape1.nelements();++i) {
1327 0 : Int minLength = shape1[i];
1328 0 : if (shape2[i]<minLength) {
1329 0 : minLength = shape2[i];
1330 : }
1331 0 : AlwaysAssert(minLength>=0, AipsError);
1332 : //if (minLength % 2 != 0) {
1333 : // if the number of pixels is odd, ensure that the centre stays
1334 : // the same by making this number even
1335 : //--minLength; // this code is a mistake and should be removed
1336 : //}
1337 0 : const Int increment1 = shape1[i] - minLength;
1338 0 : const Int increment2 = shape2[i] - minLength;
1339 0 : blc1[i] += increment1/2;
1340 0 : trc1[i] -= increment1/2 + (increment1 % 2 != 0 ? 1 : 0);
1341 0 : blc2[i] += increment2/2;
1342 0 : trc2[i] -= increment2/2 + (increment2 % 2 != 0 ? 1 : 0);
1343 : }
1344 0 : }
1345 :
1346 0 : Int MatrixCleaner::findBeamPatch(const Float maxScaleSize, const Int& nx, const Int& ny,
1347 : const Float psfBeam, const Float nBeams)
1348 : {
1349 0 : Int psupport = (Int) ( sqrt( psfBeam*psfBeam + maxScaleSize*maxScaleSize ) * nBeams );
1350 :
1351 : // At least this big...
1352 0 : if(psupport < psfBeam*nBeams ) psupport = static_cast<Int>(psfBeam*nBeams);
1353 :
1354 : // Not too big...
1355 0 : if(psupport > nx || psupport > ny) psupport = std::min(nx,ny);
1356 :
1357 : // Make it even.
1358 0 : if (psupport%2 != 0) psupport -= 1;
1359 :
1360 0 : return psupport;
1361 : }
1362 :
1363 :
1364 : } //# NAMESPACE CASA - END
|
