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: aips2-request@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 215 : Bool MatrixCleaner::validatePsf(const Matrix<Float> & psf)
67 : {
68 430 : LogIO os(LogOrigin("MatrixCleaner", "validatePsf()", WHERE));
69 :
70 : // Find the peak of the raw Psf
71 215 : AlwaysAssert(psf.shape().product() != 0, AipsError);
72 215 : Float maxPsf=0;
73 215 : itsPositionPeakPsf=IPosition(psf.shape().nelements(), 0);
74 : //findMaxAbs(psf, maxPsf, itsPositionPeakPsf);
75 215 : Int psfSupport = findBeamPatch(0.0, psf.shape()(0), psf.shape()(1),
76 215 : 4.0, 20.0);
77 215 : findPSFMaxAbs(psf, maxPsf, itsPositionPeakPsf, psfSupport);
78 215 : os << "Peak of PSF = " << maxPsf << " at " << itsPositionPeakPsf
79 215 : << LogIO::POST;
80 430 : return true;
81 : }
82 :
83 :
84 187 : MatrixCleaner::MatrixCleaner():
85 : itsMask( ),
86 : itsSmallScaleBias(0.0),
87 : itsMaskThreshold(0.9),
88 : itsDirty( ),
89 : itsXfr( ),
90 : itsScaleSizes(0),
91 : itsMaximumResidual(0.0),
92 : itsStrengthOptimum(0.0),
93 : itsTotalFlux(0.0),
94 : itsChoose(true),
95 : itsDoSpeedup(false),
96 : itsIgnoreCenterBox(false),
97 : itsStopAtLargeScaleNegative(false),
98 : itsStopPointMode(-1),
99 : itsDidStopPointMode(false),
100 : itsJustStarting(true),
101 : psfShape_p(0),
102 187 : noClean_p(false)
103 : {
104 187 : itsMemoryMB=Double(HostInfo::memoryTotal()/1024)/16.0;
105 187 : itsScales.resize(0);
106 187 : itsScaleXfrs.resize(0);
107 187 : itsDirtyConvScales.resize(0);
108 187 : itsPsfConvScales.resize(0);
109 187 : itsScaleMasks.resize(0);
110 187 : itsScalesValid = false;
111 187 : itsStartingIter = 0;
112 187 : itsFracThreshold=Quantity(0.0, "%");
113 187 : }
114 :
115 :
116 :
117 0 : MatrixCleaner::MatrixCleaner(const Matrix<Float> & psf,
118 0 : const Matrix<Float> &dirty):
119 : itsMask( ),
120 : itsSmallScaleBias(0.0),
121 : itsScaleSizes(0),
122 : itsMaximumResidual(0.0),
123 : itsStrengthOptimum(0.),
124 : itsTotalFlux(0.0),
125 : itsChoose(true),
126 : itsDoSpeedup(false),
127 : itsIgnoreCenterBox(false),
128 : itsStopAtLargeScaleNegative(false),
129 : itsStopPointMode(-1),
130 : itsDidStopPointMode(false),
131 : 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 84 : void MatrixCleaner::setPsf(const Matrix<Float>& psf){
167 84 : itsXfr=new Matrix<Complex>();
168 84 : AlwaysAssert(validatePsf(psf), AipsError);
169 84 : psfShape_p.resize(0, false);
170 84 : psfShape_p=psf.shape();
171 168 : FFTServer<Float,Complex> fft(psf.shape());
172 84 : fft.fft0(*itsXfr, psf);
173 : //cout << "shapes " << itsXfr->shape() << " psf " << psf.shape() << endl;
174 84 : }
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 187 : MatrixCleaner::~MatrixCleaner()
212 : {
213 187 : destroyScales();
214 187 : if(!itsMask.null()) itsMask=0;
215 187 : }
216 :
217 :
218 44 : void MatrixCleaner::makeDirtyScales(){
219 :
220 44 : if(!itsScalesValid || itsNscales < 1 || itsDirty.null() || (itsNscales != Int(itsScaleXfrs.nelements())) )
221 0 : return;
222 :
223 44 : 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 44 : 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 88 : Matrix<Complex> dirtyFT;
233 88 : FFTServer<Float,Complex> fft(itsDirty->shape());
234 44 : fft.fft0(dirtyFT, *itsDirty);
235 44 : itsDirtyConvScales.resize(itsNscales, true);
236 44 : 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 256 : for (scale=0; scale<itsNscales; scale++) {
243 424 : Matrix<Complex> cWork;
244 : // Dirty * scale
245 : // cout << "scale " << scale << " itsScaleptr " << &(itsScaleXfrs[scale]) << "\n"<< endl;
246 :
247 212 : itsDirtyConvScales[scale]=Matrix<Float>(itsDirty->shape());
248 212 : cWork=((dirtyFT)*(itsScaleXfrs[scale]));
249 212 : fft.fft0((itsDirtyConvScales[scale]), cWork, false);
250 212 : fft.flip((itsDirtyConvScales[scale]), false, false);
251 : }
252 : }
253 :
254 : }
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 44 : void MatrixCleaner::setMask(Matrix<Float> & mask, const Float& maskThreshold)
271 : {
272 44 : itsMaskThreshold = maskThreshold;
273 88 : IPosition maskShape = mask.shape();
274 :
275 : //cerr << "Mask Shape " << mask.shape() << endl;
276 : // This is not needed after the first steps
277 44 : itsMask = new Matrix<Float>(mask.shape());
278 44 : itsMask->assign(mask);
279 44 : noClean_p=(max(*itsMask) < itsMaskThreshold) ? true : false;
280 :
281 :
282 :
283 44 : if (itsScalesValid) {
284 44 : makeScaleMasks();
285 : }
286 :
287 44 : }
288 :
289 88 : Bool MatrixCleaner::setcontrol(CleanEnums::CleanType cleanType,
290 : const Int niter,
291 : const Float gain,
292 : const Quantity& threshold)
293 : {
294 88 : return setcontrol(cleanType, niter, gain, threshold, Quantity(0.0, "%"));
295 : }
296 :
297 : // Set up the control parameters
298 88 : Bool MatrixCleaner::setcontrol(CleanEnums::CleanType cleanType,
299 : const Int niter,
300 : const Float gain,
301 : const Quantity& aThreshold,
302 : const Quantity& fThreshold)
303 : {
304 88 : itsCleanType=cleanType;
305 88 : itsMaxNiter=niter;
306 88 : itsGain=gain;
307 88 : itsThreshold=aThreshold;
308 88 : itsFracThreshold=fThreshold;
309 88 : 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 44 : Int MatrixCleaner::clean(Matrix<Float>& model,
323 : Bool /*showProgress*/)
324 : {
325 44 : AlwaysAssert(model.shape()==itsDirty->shape(), AipsError);
326 :
327 132 : LogIO os(LogOrigin("MatrixCleaner", "clean()", WHERE));
328 :
329 44 : Float tmpMaximumResidual=0.0;
330 :
331 44 : Int nScalesToClean=itsNscales;
332 44 : if (itsCleanType==CleanEnums::HOGBOM) {
333 0 : os << LogIO::NORMAL1 << "Hogbom clean algorithm" << LogIO::POST;
334 0 : nScalesToClean=1;
335 : }
336 44 : else if (itsCleanType==CleanEnums::MULTISCALE) {
337 44 : if (nScalesToClean==1) {
338 9 : os << LogIO::NORMAL1 << "Multi-scale clean with only one scale" << LogIO::POST;
339 : }
340 : else {
341 35 : os << LogIO::NORMAL1 << "Multi-scale clean algorithm" << LogIO::POST;
342 : }
343 : }
344 :
345 : Int scale;
346 88 : Vector<Float> scaleBias(nScalesToClean);
347 44 : if (nScalesToClean > 1) {
348 238 : for (scale=0;scale<nScalesToClean;scale++) {
349 812 : scaleBias(scale) = 1 - itsSmallScaleBias *
350 203 : itsScaleSizes(scale)/itsScaleSizes(nScalesToClean-1);
351 203 : os << "scale " << scale+1 << " = " << itsScaleSizes(scale) << " pixels with bias = " << scaleBias(scale) << LogIO::POST;
352 : }
353 : } else {
354 9 : scaleBias(0) = 1.0;
355 : }
356 :
357 44 : AlwaysAssert(itsScalesValid, AipsError);
358 : ////no need to use all cores if possible
359 44 : Int nth=nScalesToClean;
360 : #ifdef _OPENMP
361 :
362 44 : nth=min(nth, omp_get_max_threads());
363 :
364 : #endif
365 : // Find the peaks of the convolved Psfs
366 88 : Vector<Float> maxPsfConvScales(nScalesToClean);
367 44 : Int naxes=model.shape().nelements();
368 44 : #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 253 : for (scale=0;scale<nScalesToClean;scale++) {
383 211 : if ( maxPsfConvScales(scale) < 0.0) {
384 : os << "As Peak of PSF is negative, you should setscales again with a smaller scale size"
385 2 : << LogIO::SEVERE;
386 2 : return -1;
387 : }
388 : }
389 :
390 : // Define a subregion for the inner quarter
391 84 : IPosition blcDirty(model.shape().nelements(), 0);
392 84 : IPosition trcDirty(model.shape()-1);
393 :
394 42 : if(!itsMask.null()){
395 42 : os << "Cleaning using given mask" << LogIO::POST;
396 42 : 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 42 : if(noClean_p)
403 0 : return 0;
404 : os << LogIO::NORMAL
405 42 : << "Cleaning pixels with mask values above " << itsMaskThreshold
406 42 : << LogIO::POST;
407 : }
408 :
409 :
410 42 : Int nx=model.shape()(0);
411 42 : Int ny=model.shape()(1);
412 :
413 :
414 42 : AlwaysAssert(itsMask->shape()(0)==nx, AipsError);
415 42 : AlwaysAssert(itsMask->shape()(1)==ny, AipsError);
416 42 : Int xbeg=nx-1;
417 42 : Int ybeg=ny-1;
418 42 : Int xend=0;
419 42 : Int yend=0;
420 23602 : for (Int iy=0;iy<ny;iy++) {
421 25074560 : for (Int ix=0;ix<nx;ix++) {
422 25051000 : if((*itsMask)(ix,iy)>0.000001) {
423 1446960 : xbeg=min(xbeg,ix);
424 1446960 : ybeg=min(ybeg,iy);
425 1446960 : xend=max(xend,ix);
426 1446960 : yend=max(yend,iy);
427 : }
428 : }
429 : }
430 :
431 42 : 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 42 : blcDirty(0)=xbeg;
447 42 : blcDirty(1)=ybeg;
448 42 : trcDirty(0)=xend;
449 42 : 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 84 : LCBox centerBox(blcDirty, trcDirty, model.shape());
466 :
467 84 : Block<Matrix<Float> > scaleMaskSubs;
468 42 : if (!itsMask.null()) {
469 42 : scaleMaskSubs.resize(itsNscales);
470 240 : for (Int is=0; is < itsNscales; is++) {
471 198 : scaleMaskSubs[is] = ((itsScaleMasks[is]))(blcDirty, trcDirty);
472 : }
473 : }
474 :
475 : // Start the iteration
476 84 : Vector<Float> maxima(nScalesToClean);
477 84 : Block<IPosition> posMaximum(nScalesToClean);
478 84 : Vector<Float> totalFluxScale(nScalesToClean); totalFluxScale=0.0;
479 42 : Float totalFlux=0.0;
480 42 : Int converged=0;
481 42 : Int stopPointModeCounter = 0;
482 42 : Int optimumScale=0;
483 42 : itsStrengthOptimum=0.0;
484 84 : IPosition positionOptimum(model.shape().nelements(), 0);
485 42 : os << "Starting iteration"<< LogIO::POST;
486 :
487 : //
488 42 : Int nx=model.shape()(0);
489 42 : Int ny=model.shape()(1);
490 84 : IPosition gip;
491 42 : gip = IPosition(2,nx,ny);
492 42 : casacore::Block<casacore::Matrix<casacore::Float> > vecWork_p;
493 42 : vecWork_p.resize(nScalesToClean);
494 :
495 240 : for(Int i=0;i<nScalesToClean;i++)
496 : {
497 198 : vecWork_p[i].resize(gip);
498 : }
499 : //
500 :
501 42 : itsIteration = itsStartingIter;
502 1945 : for (Int ii=itsStartingIter; ii < itsMaxNiter; ii++) {
503 1910 : itsIteration++;
504 : // Find the peak residual
505 1910 : itsStrengthOptimum = 0.0;
506 1910 : optimumScale = 0;
507 :
508 1910 : #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 12150 : for (scale=0; scale<nScalesToClean; scale++) {
540 10240 : if(abs(maxima(scale))>abs(itsStrengthOptimum)) {
541 4491 : optimumScale=scale;
542 4491 : itsStrengthOptimum=maxima(scale);
543 4491 : 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 1910 : 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 1910 : 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 1910 : itsStrengthOptimum /= scaleBias(optimumScale);
555 1910 : itsStrengthOptimum /= (itsDirtyConvScales[optimumScale])(posMaximum[optimumScale]);
556 : }
557 :
558 1910 : AlwaysAssert(optimumScale<nScalesToClean, AipsError);
559 :
560 : // Now add to the total flux
561 1910 : totalFlux += (itsStrengthOptimum*itsGain);
562 1910 : itsTotalFlux=totalFlux;
563 1910 : totalFluxScale(optimumScale) += (itsStrengthOptimum*itsGain);
564 :
565 1910 : if(ii==itsStartingIter ) {
566 42 : itsMaximumResidual=abs(itsStrengthOptimum);
567 42 : tmpMaximumResidual=itsMaximumResidual;
568 42 : os << "Initial maximum residual is " << itsMaximumResidual;
569 42 : if( !itsMask.null() ) { os << " within the mask "; }
570 42 : os << LogIO::POST;
571 : }
572 :
573 : // Various ways of stopping:
574 : // 1. stop if below threshold
575 1910 : if(abs(itsStrengthOptimum)<threshold() ) {
576 8 : os << "Reached stopping threshold " << threshold() << " at iteration "<<
577 4 : ii << LogIO::POST;
578 4 : os << "Optimum flux is " << abs(itsStrengthOptimum) << LogIO::POST;
579 4 : converged = 1;
580 7 : break;
581 : }
582 : // 2. negatives on largest scale?
583 1906 : 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 1906 : 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 1906 : if(((abs(itsStrengthOptimum)-abs(tmpMaximumResidual)) > (abs(tmpMaximumResidual)/2.0))
609 1906 : && !(itsStopAtLargeScaleNegative)){
610 : os << "Diverging due to large scale?"
611 3 : << LogIO::POST;
612 : //clean is diverging most probably due to the large scale
613 3 : converged=-2;
614 3 : break;
615 : }
616 : //5. Diverging for some other reason; may just need another CS-style reconciling
617 1903 : 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 1903 : if (itsIteration == itsStartingIter + 1) {
635 42 : os << "iteration MaximumResidual CleanedFlux" << LogIO::POST;
636 : }
637 1903 : if ((itsIteration % (itsMaxNiter/10 > 0 ? itsMaxNiter/10 : 1)) == 0) {
638 : //Good place to re-up the fiducial maximum residual
639 : //tmpMaximumResidual=abs(itsStrengthOptimum);
640 331 : os << itsIteration <<" "<<itsStrengthOptimum<<" "
641 331 : << totalFlux <<LogIO::POST;
642 : }
643 : }
644 :
645 : Float scaleFactor;
646 1903 : 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 3806 : IPosition support(model.shape());
651 1903 : support(0)=max(Int(itsScaleSizes(itsNscales-1)+0.5), support(0));
652 1903 : support(1)=max(Int(itsScaleSizes(itsNscales-1)+0.5), support(1));
653 :
654 3806 : 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 3806 : IPosition blc(positionOptimum-support/2);
666 5709 : IPosition trc(positionOptimum+support/2-1);
667 1903 : LCBox::verify(blc, trc, inc, model.shape());
668 :
669 : //cout << "blc " << blc.asVector() << " trc " << trc.asVector() << endl;
670 :
671 3806 : IPosition blcPsf(blc+itsPositionPeakPsf-positionOptimum);
672 3806 : IPosition trcPsf(trc+itsPositionPeakPsf-positionOptimum);
673 1903 : LCBox::verify(blcPsf, trcPsf, inc, model.shape());
674 1903 : 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 3806 : Matrix<Float> modelSub=model(blc, trc);
681 3806 : Matrix<Float> scaleSub=(itsScales[optimumScale])(blcPsf,trcPsf);
682 :
683 :
684 : // Now do the addition of this scale to the model image....
685 1903 : modelSub += scaleFactor*scaleSub;
686 :
687 1903 : #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 1903 : blcDirty = blc;
701 1903 : trcDirty = trc;
702 : }
703 : // End of iteration
704 :
705 240 : for (scale=0;scale<nScalesToClean;scale++) {
706 : os << LogIO::NORMAL
707 396 : << " " << scale << " " << totalFluxScale(scale)
708 198 : << 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 42 : if(!converged) {
721 35 : os << "Failed to reach stopping threshold" << LogIO::POST;
722 : }
723 :
724 42 : return converged;
725 : }
726 :
727 215 : Bool MatrixCleaner::findPSFMaxAbs(const Matrix<Float>& lattice,
728 : Float& maxAbs,
729 : IPosition& posMaxAbs,
730 : const Int& supportSize)
731 : {
732 645 : LogIO os(LogOrigin("MatrixCleaner", "findPSFMaxAbs()", WHERE));
733 215 : posMaxAbs = IPosition(lattice.shape().nelements(), 0);
734 215 : maxAbs=0.0;
735 :
736 215 : Float maxVal=0.0;
737 215 : IPosition psf2DShape(lattice.shape());
738 215 : Int blc0 = (psf2DShape(0) > supportSize) ? psf2DShape(0)/2 - supportSize/2 : 0;
739 215 : Int blc1 = (psf2DShape(1) > supportSize) ? psf2DShape(1)/2 - supportSize/2 : 0;
740 215 : Int trc0 = (psf2DShape(0) > supportSize) ? psf2DShape(0)/2 + supportSize/2 : psf2DShape(0)-1;
741 :
742 215 : 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 15869 : for (Int j=blc1; j < trc1; ++j)
748 1252560 : for (Int i=blc0 ; i < trc0; ++i)
749 1236906 : if ((maxAbs = abs(lattice(i,j))) > maxVal)
750 : {
751 15294 : maxVal = maxAbs;
752 15294 : posMaxAbs(0)=i; posMaxAbs(1)=j;
753 : }
754 215 : maxAbs=maxVal;
755 : //cerr << "######## " << posMaxAbs << " " << maxVal << endl;
756 430 : return true;
757 : }
758 :
759 212 : Bool MatrixCleaner::findMaxAbs(const Matrix<Float>& lattice,
760 : Float& maxAbs,
761 : IPosition& posMaxAbs)
762 : {
763 :
764 212 : posMaxAbs = IPosition(lattice.shape().nelements(), 0);
765 212 : maxAbs=0.0;
766 :
767 : Float minVal;
768 212 : IPosition posmin(lattice.shape().nelements(), 0);
769 212 : minMax(minVal, maxAbs, posmin, posMaxAbs, lattice);
770 : //cout << "min " << minVal << " " << maxAbs << " " << max(lattice) << endl;
771 212 : if(abs(minVal) > abs(maxAbs)){
772 3 : maxAbs=minVal;
773 3 : posMaxAbs=posmin;
774 : }
775 424 : return true;
776 : }
777 :
778 :
779 :
780 :
781 :
782 25944 : Bool MatrixCleaner::findMaxAbsMask(const Matrix<Float>& lattice,
783 : const Matrix<Float>& mask,
784 : Float& maxAbs,
785 : IPosition& posMaxAbs)
786 : {
787 :
788 25944 : posMaxAbs = IPosition(lattice.shape().nelements(), 0);
789 25944 : maxAbs=0.0;
790 : Float minVal;
791 25944 : IPosition posmin(lattice.shape().nelements(), 0);
792 25944 : minMaxMasked(minVal, maxAbs, posmin, posMaxAbs, lattice, mask);
793 25944 : if(abs(minVal) > abs(maxAbs)){
794 2253 : maxAbs=minVal;
795 2253 : posMaxAbs=posmin;
796 : }
797 :
798 51888 : return true;
799 : }
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 : }
830 :
831 :
832 84 : void MatrixCleaner::setDirty(const Matrix<Float>& dirty){
833 :
834 84 : itsDirty=new Matrix<Float>(dirty.shape());
835 84 : itsDirty->assign(dirty);
836 :
837 :
838 :
839 84 : }
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 44 : void MatrixCleaner::defineScales(const Vector<Float>& scaleSizes){
845 44 : if(itsScales.nelements()>0) {
846 21 : destroyScales();
847 : }
848 :
849 44 : destroyMasks();
850 44 : itsNscales=scaleSizes.nelements();
851 44 : itsScaleSizes.resize(itsNscales);
852 44 : itsScaleSizes=scaleSizes; // make a copy that we can call our own
853 44 : GenSort<Float>::sort(itsScaleSizes);
854 44 : itsScalesValid=false;
855 44 : }
856 :
857 44 : void MatrixCleaner::makePsfScales(){
858 132 : LogIO os(LogOrigin("MatrixCleaner", "makePsfScales()", WHERE));
859 44 : if(itsNscales < 1)
860 0 : throw(AipsError("Scales have to be set"));
861 44 : if(itsXfr.null())
862 0 : throw(AipsError("Psf is not defined"));
863 44 : destroyScales();
864 44 : itsScales.resize(itsNscales, true);
865 44 : itsScaleXfrs.resize(itsNscales, true);
866 44 : itsPsfConvScales.resize((itsNscales+1)*(itsNscales+1), true);
867 88 : FFTServer<Float,Complex> fft(psfShape_p);
868 44 : Int scale=0;
869 256 : for(scale=0; scale<itsNscales;scale++) {
870 212 : itsScales[scale] = Matrix<Float>(psfShape_p);
871 212 : makeScale(itsScales[scale], itsScaleSizes(scale));
872 212 : itsScaleXfrs[scale] = Matrix<Complex> ();
873 212 : fft.fft0(itsScaleXfrs[scale], itsScales[scale]);
874 : }
875 44 : Matrix<Complex> cWork;
876 :
877 256 : for (scale=0; scale<itsNscales;scale++) {
878 212 : os << "Calculating convolutions for scale " << scale << LogIO::POST;
879 : //PSF * scale
880 212 : itsPsfConvScales[scale] = Matrix<Float>(psfShape_p);
881 212 : cWork=((*itsXfr)*(itsScaleXfrs[scale])*(itsScaleXfrs[scale]));
882 : //cout << "shape " << cWork.shape() << " " << itsPsfConvScales[scale].shape() << endl;
883 :
884 212 : fft.fft0((itsPsfConvScales[scale]), cWork, false);
885 212 : fft.flip(itsPsfConvScales[scale], false, false);
886 :
887 : //cout << "psf scale " << scale << " " << max(itsPsfConvScales[scale]) << " " << min(itsPsfConvScales[scale]) << endl;
888 :
889 949 : for (Int otherscale=scale;otherscale<itsNscales;otherscale++) {
890 :
891 737 : AlwaysAssert(index(scale, otherscale)<Int(itsPsfConvScales.nelements()),
892 : AipsError);
893 :
894 : // PSF * scale * otherscale
895 737 : itsPsfConvScales[index(scale,otherscale)] =Matrix<Float>(psfShape_p);
896 737 : cWork=((*itsXfr)*(itsScaleXfrs[scale])*(itsScaleXfrs[otherscale]));
897 737 : 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 44 : itsScalesValid=true;
904 :
905 44 : }
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 : }
1020 :
1021 : // Make a single scale size image
1022 416 : void MatrixCleaner::makeScale(Matrix<Float>& iscale, const Float& scaleSize)
1023 : {
1024 :
1025 416 : Int nx=iscale.shape()(0);
1026 416 : Int ny=iscale.shape()(1);
1027 : //Matrix<Float> iscale(nx, ny);
1028 416 : iscale=0.0;
1029 :
1030 416 : Double refi=nx/2;
1031 416 : Double refj=ny/2;
1032 :
1033 416 : if(scaleSize==0.0) {
1034 163 : iscale(Int(refi), Int(refj)) = 1.0;
1035 : }
1036 : else {
1037 253 : AlwaysAssert(scaleSize>0.0,AipsError);
1038 :
1039 253 : Int mini = max( 0, (Int)(refi-scaleSize));
1040 253 : Int maxi = min(nx-1, (Int)(refi+scaleSize));
1041 253 : Int minj = max( 0, (Int)(refj-scaleSize));
1042 253 : Int maxj = min(ny-1, (Int)(refj+scaleSize));
1043 :
1044 253 : Float ypart=0.0;
1045 253 : Float volume=0.0;
1046 253 : Float rad2=0.0;
1047 253 : Float rad=0.0;
1048 :
1049 22814 : for (Int j=minj;j<=maxj;j++) {
1050 22561 : ypart = square( (refj - (Double)(j)) / scaleSize );
1051 3010828 : for (Int i=mini;i<=maxi;i++) {
1052 2988267 : rad2 = ypart + square( (refi - (Double)(i)) / scaleSize );
1053 2988267 : if (rad2 < 1.0) {
1054 2315494 : if (rad2 <= 0.0) {
1055 253 : rad = 0.0;
1056 : } else {
1057 2315241 : rad = sqrt(rad2);
1058 : }
1059 2315494 : iscale(i,j) = (1.0 - rad2) * spheroidal(rad);
1060 2315494 : volume += iscale(i,j);
1061 : } else {
1062 672773 : iscale(i,j) = 0.0;
1063 : }
1064 : }
1065 : }
1066 253 : iscale/=volume;
1067 : }
1068 : //scale.putSlice(iscale, IPosition(scale.ndim(),0), IPosition(scale.ndim(),1));
1069 416 : }
1070 :
1071 : // Calculate the spheroidal function
1072 2315494 : Float MatrixCleaner::spheroidal(Float nu) {
1073 :
1074 2315494 : if (nu <= 0) {
1075 253 : return 1.0;
1076 2315241 : } else if (nu >= 1.0) {
1077 0 : return 0.0;
1078 : } else {
1079 2315241 : uInt np = 5;
1080 2315241 : uInt nq = 3;
1081 4630482 : Matrix<Float> p(np, 2);
1082 4630482 : Matrix<Float> q(nq, 2);
1083 2315241 : p(0,0) = 8.203343e-2;
1084 2315241 : p(1,0) = -3.644705e-1;
1085 2315241 : p(2,0) = 6.278660e-1;
1086 2315241 : p(3,0) = -5.335581e-1;
1087 2315241 : p(4,0) = 2.312756e-1;
1088 2315241 : p(0,1) = 4.028559e-3;
1089 2315241 : p(1,1) = -3.697768e-2;
1090 2315241 : p(2,1) = 1.021332e-1;
1091 2315241 : p(3,1) = -1.201436e-1;
1092 2315241 : p(4,1) = 6.412774e-2;
1093 2315241 : q(0,0) = 1.0000000e0;
1094 2315241 : q(1,0) = 8.212018e-1;
1095 2315241 : q(2,0) = 2.078043e-1;
1096 2315241 : q(0,1) = 1.0000000e0;
1097 2315241 : q(1,1) = 9.599102e-1;
1098 2315241 : q(2,1) = 2.918724e-1;
1099 2315241 : uInt part = 0;
1100 2315241 : Float nuend = 0.0;
1101 2315241 : if (nu >= 0.0 && nu < 0.75) {
1102 1314217 : part = 0;
1103 1314217 : nuend = 0.75;
1104 1001024 : } else if (nu >= 0.75 && nu <= 1.00) {
1105 1001024 : part = 1;
1106 1001024 : nuend = 1.0;
1107 : }
1108 :
1109 2315241 : Float top = p(0,part);
1110 2315241 : Float delnusq = pow(nu,2.0) - pow(nuend,2.0);
1111 : uInt k;
1112 11576205 : for (k=1; k<np; k++) {
1113 9260964 : top += p(k, part) * pow(delnusq, (Float)k);
1114 : }
1115 2315241 : Float bot = q(0, part);
1116 6945723 : for (k=1; k<nq; k++) {
1117 4630482 : bot += q(k,part) * pow(delnusq, (Float)k);
1118 : }
1119 :
1120 2315241 : if (bot != 0.0) {
1121 2315241 : return (top/bot);
1122 : } else {
1123 0 : return 0.0;
1124 : }
1125 : }
1126 : }
1127 :
1128 :
1129 : // Calculate index into PsfConvScales
1130 12418 : Int MatrixCleaner::index(const Int scale, const Int otherscale) {
1131 12418 : if(otherscale>scale) {
1132 5851 : return scale + itsNscales*(otherscale+1);
1133 : }
1134 : else {
1135 6567 : return otherscale + itsNscales*(scale+1);
1136 : }
1137 : }
1138 :
1139 :
1140 288 : Bool MatrixCleaner::destroyScales()
1141 : {
1142 288 : if(!itsScalesValid) return true;
1143 288 : for(uInt scale=0; scale<itsScales.nelements();scale++) {
1144 : //if(itsScales[scale]) delete itsScales[scale];
1145 212 : itsScales[scale].resize();
1146 : }
1147 288 : for(uInt scale=0; scale<itsScaleXfrs.nelements();scale++) {
1148 : //if(itsScaleXfrs[scale]) delete itsScaleXfrs[scale];
1149 212 : itsScaleXfrs[scale].resize();
1150 : }
1151 288 : for(uInt scale=0; scale<itsDirtyConvScales.nelements();scale++) {
1152 : //if(itsDirtyConvScales[scale]) delete itsDirtyConvScales[scale];
1153 212 : itsDirtyConvScales[scale].resize();
1154 : }
1155 1806 : for(uInt scale=0; scale<itsPsfConvScales.nelements();scale++) {
1156 : //if(itsPsfConvScales[scale]) delete itsPsfConvScales[scale];
1157 1730 : itsPsfConvScales[scale].resize();
1158 : }
1159 76 : destroyMasks();
1160 76 : itsScales.resize(0, true);
1161 76 : itsDirtyConvScales.resize(0,true);
1162 76 : itsPsfConvScales.resize(0, true);
1163 76 : itsScalesValid=false;
1164 76 : return true;
1165 : }
1166 :
1167 :
1168 :
1169 204 : Bool MatrixCleaner::destroyMasks()
1170 : {
1171 646 : for(uInt scale=0; scale<itsScaleMasks.nelements();scale++) {
1172 : //if(itsScaleMasks[scale]) delete itsScaleMasks[scale];
1173 442 : itsScaleMasks[scale].resize();
1174 : }
1175 204 : itsScaleMasks.resize(0);
1176 204 : 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 44 : Bool MatrixCleaner::makeScaleMasks()
1195 : {
1196 132 : LogIO os(LogOrigin("MatrixCleaner", "makeScaleMasks()", WHERE));
1197 : Int scale;
1198 :
1199 :
1200 44 : if(!itsScalesValid) {
1201 : os << "Scales are not yet set - cannot set scale masks"
1202 0 : << LogIO::EXCEPTION;
1203 : }
1204 :
1205 :
1206 44 : destroyMasks();
1207 :
1208 44 : if(itsMask.null() || noClean_p)
1209 0 : return false;
1210 :
1211 44 : AlwaysAssert((itsMask->shape() == psfShape_p), AipsError);
1212 : //cout << "itsmask " << itsMask->shape() << endl;
1213 88 : FFTServer<Float,Complex> fft(itsMask->shape());
1214 :
1215 88 : Matrix<Complex> maskFT;
1216 44 : fft.fft0(maskFT, *itsMask);
1217 44 : itsScaleMasks.resize(itsNscales, true);
1218 : // Now we can do all the convolutions
1219 44 : Matrix<Complex> cWork;
1220 256 : for (scale=0; scale<itsNscales;scale++) {
1221 : //AlwaysAssert(itsScaleXfrs[scale], AipsError);
1222 212 : 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 212 : itsScaleMasks[scale] = Matrix<Float>(itsMask->shape());
1228 212 : cWork=((maskFT)*(itsScaleXfrs[scale]));
1229 212 : fft.fft0(itsScaleMasks[scale], cWork, false);
1230 212 : fft.flip(itsScaleMasks[scale], false, false);
1231 167652 : for (Int j=0 ; j < (itsMask->shape())(1); ++j){
1232 193280040 : for (Int k =0 ; k < (itsMask->shape())(0); ++k){
1233 193112600 : if(itsMaskThreshold > 0){
1234 193112600 : (itsScaleMasks[scale])(k,j) = (itsScaleMasks[scale])(k,j) > itsMaskThreshold ? 1.0 : 0.0;
1235 : }
1236 : }
1237 : }
1238 212 : Float mysum = sum(itsScaleMasks[scale] );
1239 212 : 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 44 : Int nx=itsScaleMasks[0].shape()(0);
1247 44 : 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 256 : for(Int scale=0;scale<itsNscales;scale++)
1252 : {
1253 212 : Int border = (Int)(itsScaleSizes[scale]*1.5);
1254 : // bottom
1255 424 : IPosition blc1(2, 0 , 0 );
1256 424 : IPosition trc1(2,nx-1, border );
1257 212 : IPosition inc1(2, 1);
1258 212 : LCBox::verify(blc1,trc1,inc1,itsScaleMasks[scale].shape());
1259 212 : (itsScaleMasks[scale])(blc1,trc1) = 0.0;
1260 : // top
1261 212 : blc1[0]=0; blc1[1]=ny-border-1;
1262 212 : trc1[0]=nx-1; trc1[1]=ny-1;
1263 212 : LCBox::verify(blc1,trc1,inc1,itsScaleMasks[scale].shape());
1264 212 : (itsScaleMasks[scale])(blc1,trc1) = 0.0;
1265 : // left
1266 212 : blc1[0]=0; blc1[1]=border;
1267 212 : trc1[0]=border; trc1[1]=ny-border-1;
1268 212 : LCBox::verify(blc1,trc1,inc1,itsScaleMasks[scale].shape());
1269 212 : (itsScaleMasks[scale])(blc1,trc1) = 0.0;
1270 : // right
1271 212 : blc1[0]=nx-border-1; blc1[1]=border;
1272 212 : trc1[0]=nx; trc1[1]=ny-border-1;
1273 212 : LCBox::verify(blc1,trc1,inc1,itsScaleMasks[scale].shape());
1274 212 : (itsScaleMasks[scale])(blc1,trc1) = 0.0;
1275 : }
1276 :
1277 44 : return true;
1278 : }
1279 :
1280 :
1281 44 : Matrix<casacore::Float> MatrixCleaner::residual(const Matrix<Float>& model){
1282 88 : FFTServer<Float,Complex> fft(model.shape());
1283 44 : if(!itsDirty.null() && (model.shape() != (itsDirty->shape())))
1284 0 : throw(AipsError("MatrixCleaner: model size has to be consistent with dirty image size"));
1285 44 : if(itsXfr.null())
1286 0 : throw(AipsError("MatrixCleaner: need to know the psf to calculate the residual"));
1287 :
1288 88 : Matrix<Complex> modelFT;
1289 : ///Convolve model with psf
1290 44 : fft.fft0(modelFT, model);
1291 44 : modelFT= (*itsXfr) * modelFT;
1292 44 : Matrix<Float> newResidual(itsDirty->shape());
1293 44 : fft.fft0(newResidual, modelFT, false);
1294 44 : fft.flip(newResidual, false, false);
1295 44 : newResidual=(*itsDirty)-newResidual;
1296 88 : return newResidual;
1297 : }
1298 :
1299 :
1300 3047 : Float MatrixCleaner::threshold() const
1301 : {
1302 3047 : if (! itsDoSpeedup) {
1303 6094 : return max(itsFracThreshold.get("%").getValue() * itsMaximumResidual /100.0,
1304 9141 : 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 4640 : void MatrixCleaner::makeBoxesSameSize(IPosition& blc1, IPosition& trc1,
1316 : IPosition &blc2, IPosition& trc2)
1317 : {
1318 4640 : const IPosition shape1 = trc1 - blc1;
1319 4640 : const IPosition shape2 = trc2 - blc2;
1320 :
1321 4640 : AlwaysAssert(shape1.nelements() == shape2.nelements(), AipsError);
1322 :
1323 4640 : if (shape1 == shape2) {
1324 4640 : 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 : }
1345 :
1346 346 : Int MatrixCleaner::findBeamPatch(const Float maxScaleSize, const Int& nx, const Int& ny,
1347 : const Float psfBeam, const Float nBeams)
1348 : {
1349 346 : Int psupport = (Int) ( sqrt( psfBeam*psfBeam + maxScaleSize*maxScaleSize ) * nBeams );
1350 :
1351 : // At least this big...
1352 346 : if(psupport < psfBeam*nBeams ) psupport = static_cast<Int>(psfBeam*nBeams);
1353 :
1354 : // Not too big...
1355 346 : if(psupport > nx || psupport > ny) psupport = std::min(nx,ny);
1356 :
1357 : // Make it even.
1358 346 : if (psupport%2 != 0) psupport -= 1;
1359 :
1360 346 : return psupport;
1361 : }
1362 :
1363 :
1364 : } //# NAMESPACE CASA - END
|
