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 : #include <synthesis/MeasurementEquations/MatrixCleaner.h>
27 : #include <synthesis/MeasurementEquations/ImageMSCleaner.h>
28 : #include <casacore/images/Images/PagedImage.h>
29 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
30 : #include <casacore/coordinates/Coordinates/CoordinateUtil.h>
31 :
32 : using namespace casacore;
33 : namespace casa { //# NAMESPACE CASA - BEGIN
34 :
35 0 : ImageMSCleaner::ImageMSCleaner(): psf_p(0), dirty_p(0), mask_p(0), nPsfChan_p(0),
36 : nImChan_p(0), nPsfPol_p(0), nImPol_p(0), chanAxis_p(-1),
37 0 : polAxis_p(-1), nMaskChan_p(0), nMaskPol_p(0), maskThresh_p(0.9), maxResidual_p(0.0)
38 : {
39 :
40 :
41 0 : }
42 0 : ImageMSCleaner::ImageMSCleaner(ImageInterface<Float>& psf,
43 0 : ImageInterface<Float>& dirty): psf_p(&psf),
44 : dirty_p(&dirty), mask_p(0),
45 : nMaskChan_p(0), nMaskPol_p(0),
46 0 : maskThresh_p(0.9), maxResidual_p(0.0){
47 :
48 0 : chanAxis_p=CoordinateUtil::findSpectralAxis(dirty_p->coordinates());
49 0 : Vector<Stokes::StokesTypes> whichPols;
50 0 : polAxis_p=CoordinateUtil::findStokesAxis(whichPols, dirty_p->coordinates());
51 0 : if(chanAxis_p > -1){
52 0 : nPsfChan_p=psf_p->shape()(chanAxis_p);
53 0 : nImChan_p=dirty_p->shape()(chanAxis_p);
54 : }
55 0 : if(polAxis_p > -1){
56 0 : nPsfPol_p=psf_p->shape()(polAxis_p);
57 0 : nImPol_p=dirty_p->shape()(polAxis_p);
58 : }
59 0 : }
60 :
61 :
62 0 : ImageMSCleaner::ImageMSCleaner(const ImageMSCleaner& other){
63 0 : operator=(other);
64 0 : }
65 :
66 0 : ImageMSCleaner::~ImageMSCleaner(){
67 :
68 0 : }
69 :
70 0 : ImageMSCleaner& ImageMSCleaner::operator=(const ImageMSCleaner& other){
71 0 : if (this != &other) {
72 0 : matClean_p=other.matClean_p;
73 0 : psf_p=other.psf_p;
74 0 : dirty_p=other.dirty_p;
75 0 : mask_p=other.mask_p;
76 0 : nPsfChan_p=other.nPsfChan_p;
77 0 : nImChan_p=other.nImChan_p;
78 0 : nPsfPol_p=other.nPsfPol_p;
79 0 : nImPol_p=other.nImPol_p;
80 0 : chanAxis_p=other.chanAxis_p;
81 0 : nMaskChan_p=other.nMaskChan_p;
82 0 : nMaskPol_p=other.nMaskPol_p;
83 0 : polAxis_p=other.polAxis_p;
84 0 : scales_p=other.scales_p;
85 0 : maskThresh_p=other.maskThresh_p;
86 0 : maxResidual_p=other.maxResidual_p;
87 : }
88 0 : return *this;
89 : }
90 0 : void ImageMSCleaner::update(ImageInterface<Float>& dirty){
91 0 : dirty_p=&dirty;
92 0 : chanAxis_p=CoordinateUtil::findSpectralAxis(dirty_p->coordinates());
93 0 : Vector<Stokes::StokesTypes> whichPols;
94 0 : polAxis_p=CoordinateUtil::findStokesAxis(whichPols, dirty_p->coordinates());
95 0 : if(chanAxis_p > -1){
96 0 : nImChan_p=dirty_p->shape()(chanAxis_p);
97 : }
98 : else
99 0 : nImChan_p=0;
100 0 : if(polAxis_p > -1){
101 0 : nImPol_p=dirty_p->shape()(polAxis_p);
102 : }
103 : else
104 0 : nImPol_p=0;
105 0 : }
106 :
107 0 : void ImageMSCleaner::setPsf(ImageInterface<Float> & psf){
108 0 : psf_p=&psf;
109 0 : Int chanAxis=CoordinateUtil::findSpectralAxis(psf_p->coordinates());
110 0 : Vector<Stokes::StokesTypes> whichPols;
111 0 : Int polAxis=CoordinateUtil::findStokesAxis(whichPols, psf_p->coordinates());
112 0 : if(chanAxis > -1){
113 0 : nPsfChan_p=psf_p->shape()(chanAxis);
114 : }
115 : else
116 0 : nPsfChan_p=0;
117 0 : if(polAxis > -1){
118 0 : nPsfPol_p=psf_p->shape()(polAxis);
119 : }
120 : else
121 0 : nPsfPol_p=0;
122 0 : }
123 :
124 0 : void ImageMSCleaner::setscales(const Int nscs, const Float scaleInc){
125 0 : LogIO os(LogOrigin("ImageMSCleaner", "setscales()", WHERE));
126 0 : Int nscales=nscs;
127 :
128 0 : if(nscales<1) {
129 0 : os << "Using default of 5 scales" << LogIO::POST;
130 0 : nscales=5;
131 : }
132 :
133 0 : scales_p.resize(nscales);
134 :
135 : // Validate scales
136 0 : os << "Creating " << nscales << " scales" << LogIO::POST;
137 0 : scales_p(0) = 0.00001 * scaleInc;
138 0 : os << "scale 0 = 0.0 arcsec" << LogIO::POST;
139 0 : for (Int scale=1; scale<nscales;scale++) {
140 0 : scales_p(scale) =
141 0 : scaleInc * pow(10.0, (Float(scale)-2.0)/2.0);
142 0 : os << "scale " << scale << " = " << scales_p(scale)
143 0 : << " arcsec" << LogIO::POST;
144 : }
145 :
146 0 : }
147 0 : void ImageMSCleaner::setscales(const Vector<Float> & scales){
148 0 : scales_p.resize(scales.nelements());
149 0 : scales_p=scales;
150 0 : }
151 :
152 0 : Bool ImageMSCleaner::setcontrol(CleanEnums::CleanType cleanType, const Int niter,
153 : const Float gain, const Quantity& aThreshold,
154 : const Quantity& fThreshold){
155 0 : return matClean_p.setcontrol(cleanType, niter, gain, aThreshold, fThreshold);
156 : }
157 0 : Bool ImageMSCleaner::setcontrol(CleanEnums::CleanType cleanType, const Int niter,
158 : const Float gain, const Quantity& threshold){
159 :
160 0 : return matClean_p.setcontrol(cleanType, niter, gain, threshold);
161 :
162 : }
163 0 : Int ImageMSCleaner::iteration() const{
164 0 : return matClean_p.iteration();
165 : }
166 0 : Int ImageMSCleaner::numberIterations() const{
167 0 : return matClean_p.numberIterations();
168 : }
169 0 : void ImageMSCleaner::startingIteration(const Int starting){
170 0 : matClean_p.startingIteration(starting);
171 0 : }
172 :
173 0 : void ImageMSCleaner::setMask(ImageInterface<Float> & mask,
174 : const Float& maskThreshold){
175 :
176 0 : maskThresh_p=maskThreshold;
177 0 : mask_p=&mask;
178 0 : Int chanAxis=CoordinateUtil::findSpectralAxis(mask_p->coordinates());
179 0 : Vector<Stokes::StokesTypes> whichPols;
180 0 : Int polAxis=CoordinateUtil::findStokesAxis(whichPols, mask_p->coordinates());
181 0 : if(chanAxis > -1){
182 0 : nMaskChan_p=mask_p->shape()(chanAxis);
183 : }
184 : else
185 0 : nMaskChan_p=0;
186 0 : if(polAxis > -1){
187 0 : nMaskPol_p=mask_p->shape()(polAxis);
188 : }
189 : else
190 0 : nMaskPol_p=0;
191 :
192 :
193 0 : }
194 :
195 0 : void ImageMSCleaner::ignoreCenterBox(Bool ign){
196 0 : matClean_p.ignoreCenterBox(ign);
197 0 : }
198 0 : void ImageMSCleaner::setSmallScaleBias(const Float x){
199 0 : matClean_p.setSmallScaleBias(x);
200 0 : }
201 0 : void ImageMSCleaner::stopAtLargeScaleNegative(){
202 0 : matClean_p.stopAtLargeScaleNegative();
203 0 : }
204 0 : void ImageMSCleaner::stopPointMode(Int nStopPointMode) {
205 0 : matClean_p.stopPointMode(nStopPointMode);
206 0 : }
207 0 : Bool ImageMSCleaner::queryStopPointMode() const{
208 0 : return matClean_p.queryStopPointMode();
209 : }
210 0 : void ImageMSCleaner::speedup(const Float Ndouble){
211 0 : matClean_p.speedup(Ndouble);
212 0 : }
213 :
214 0 : Bool ImageMSCleaner::setupMatCleaner(const String& alg, const Int niter,
215 : const Float gain, const Quantity& threshold, const Quantity& fthresh){
216 :
217 0 : LogIO os(LogOrigin("ImageMSCleaner", "setupclean()", WHERE));
218 :
219 0 : String algorithm=alg;
220 0 : algorithm.downcase();
221 0 : if((algorithm=="msclean")||(algorithm=="fullmsclean" )|| (algorithm=="multiscale")
222 0 : || (algorithm=="fullmultiscale")) {
223 : //os << "Cleaning image using multi-scale algorithm" << LogIO::POST;
224 0 : if(scales_p.nelements()==0) {
225 0 : os << LogIO::SEVERE << "Scales not yet set" << LogIO::POST;
226 0 : return false;
227 : }
228 : //matClean_p.setscales(scaleSizes_p);
229 0 : matClean_p.setcontrol(CleanEnums::MULTISCALE, niter, gain, threshold, fthresh);
230 : }
231 0 : else if (algorithm=="hogbom") {
232 0 : scales_p=Vector<Float>(1,0.0);
233 0 : matClean_p.defineScales(scales_p);
234 0 : matClean_p.setcontrol(CleanEnums::HOGBOM, niter, gain, threshold, fthresh);
235 : } else {
236 0 : os << LogIO::SEVERE << "Unknown algorithm: " << algorithm << LogIO::POST;
237 0 : return false;
238 : }
239 :
240 0 : if(algorithm=="fullmsclean" || algorithm=="fullmultiscale") {
241 0 : matClean_p.ignoreCenterBox(true);
242 : }
243 0 : return true;
244 : }
245 :
246 0 : Int ImageMSCleaner::clean(ImageInterface<Float> & modelimage, const String& algorithm,
247 : const Int niter,
248 : const Float gain, const Quantity& threshold, const Quantity& fthresh, Bool /*doPlotProgress*/){
249 :
250 :
251 :
252 :
253 :
254 0 : Int result=0;
255 : ///case of single plane mask
256 : //now may be a time to set stuff scales will be done later
257 0 : if(!setupMatCleaner(algorithm, niter, gain, threshold, fthresh))
258 0 : return false;
259 0 : matClean_p.defineScales(scales_p);
260 : //cerr << "nPol " << nMaskPol_p << " " << nPsfPol_p << " " << nImPol_p << endl;
261 : //cerr << "nChan " << nMaskChan_p << " " << nPsfChan_p << " " << nImChan_p << endl;
262 0 : if( (nMaskPol_p >1) && (nMaskPol_p != nImPol_p))
263 0 : throw(AipsError("Donot know how to deal with mask that has different pol axis as Image"));
264 0 : if(mask_p && nMaskChan_p < 2 && nMaskPol_p < 2){
265 0 : Matrix<Float> maskMat;
266 0 : Array<Float> mbuf;
267 0 : mask_p->get(mbuf, true);
268 0 : maskMat.reference(mbuf);
269 :
270 0 : matClean_p.setMask(maskMat, maskThresh_p);
271 : }
272 0 : if(!psf_p)
273 0 : throw(AipsError("No PSF defined "));
274 : ///case of single plane psf
275 0 : if(nPsfChan_p < 2){
276 0 : Matrix<Float> psfMat;
277 0 : Array<Float> mbuf;
278 0 : if(nPsfPol_p > 1){
279 : ///First stokes psf is good enough
280 0 : IPosition blc(dirty_p->ndim(),0);
281 0 : IPosition trc=(dirty_p->shape()) -1;
282 0 : trc(polAxis_p)=0;
283 0 : Slicer sl(blc, trc, Slicer::endIsLast);
284 0 : psf_p->getSlice(mbuf, sl, true);
285 : }
286 : else{
287 0 : psf_p->get(mbuf, true);
288 : }
289 0 : psfMat.reference(mbuf);
290 0 : matClean_p.setPsf(psfMat);
291 0 : matClean_p.makePsfScales();
292 : //matClean_p.makeScaleMasks();
293 :
294 : }
295 :
296 :
297 :
298 : //has cube axis
299 0 : if(chanAxis_p > -1 && nImChan_p > 0){
300 0 : Int npol=1;
301 0 : if(polAxis_p > -1 && nImPol_p >0){
302 0 : npol=nImPol_p;
303 : }
304 0 : Matrix<Float> subModel;
305 0 : Bool getModel=false;
306 0 : IPosition blc(dirty_p->ndim(), 0);
307 0 : IPosition trc=dirty_p->shape() -1;
308 0 : for (Int k=0; k < nImChan_p ; ++k){
309 0 : for (Int j=0; j < npol; ++j){
310 0 : if(npol > 1 ){
311 0 : blc(polAxis_p)=j;
312 0 : trc(polAxis_p)=j;
313 : }
314 0 : blc(chanAxis_p)=k;
315 0 : trc(chanAxis_p)=k;
316 0 : Slicer sl(blc, trc, Slicer::endIsLast);
317 0 : Matrix<Float> dirtySub ;
318 0 : Array<Float> buf;
319 0 : dirty_p->getSlice(buf, sl, true);
320 0 : dirtySub.reference(buf);
321 0 : matClean_p.setDirty(dirtySub);
322 0 : Array<Float> bufMod;
323 0 : getModel=modelimage.getSlice(bufMod, sl, true);
324 0 : subModel.reference(bufMod);
325 :
326 0 : if((nPsfChan_p>1) && (k<(nPsfChan_p-1))){
327 0 : Matrix<Float> psfSub;
328 0 : Array<Float> buf1;
329 0 : psf_p->getSlice(buf1, sl, true);
330 0 : psfSub.reference(buf1);
331 0 : matClean_p.setPsf(psfSub);
332 0 : matClean_p.makePsfScales();
333 : }
334 0 : matClean_p.makeDirtyScales();
335 0 : if(mask_p !=0){
336 0 : if((nMaskChan_p >1) && (k<(nMaskChan_p-1))){
337 0 : Matrix<Float> maskSub;
338 0 : Array<Float> buf2;
339 0 : mask_p->getSlice(buf2, sl, true);
340 0 : maskSub.reference(buf2);
341 0 : matClean_p.setMask(maskSub, maskThresh_p);
342 : //Set mask above does the makeScaleMasks as psf-scale-Xfr are valid
343 : }
344 : //use one plane mask to all planes
345 0 : else if(nMaskChan_p==1){
346 0 : if((nMaskPol_p >1) && (mask_p !=0) ){
347 0 : Matrix<Float> maskSub;
348 0 : Array<Float> buf2;
349 0 : mask_p->getSlice(buf2, sl, true);
350 0 : maskSub.reference(buf2);
351 0 : matClean_p.setMask(maskSub, maskThresh_p);
352 : }
353 0 : matClean_p.makeScaleMasks();
354 : }
355 : else{
356 0 : matClean_p.unsetMask();
357 : }
358 : }
359 0 : result=matClean_p.clean(subModel, true);
360 : //Update the private flux and residuals here
361 0 : maxResidual_p=max(maxResidual_p, matClean_p.strengthOptimum());
362 0 : if(!getModel)
363 0 : modelimage.putSlice(bufMod, sl.start());
364 :
365 : }
366 0 : }
367 : }
368 : //has no channel but has pol
369 0 : else if(polAxis_p > -1 && nImPol_p > 1){
370 : //The psf should have been set above before the if loop
371 : //for this case
372 0 : Matrix<Float> subModel;
373 0 : Bool getModel=false;
374 0 : IPosition blc(dirty_p->ndim(), 0);
375 0 : IPosition trc=dirty_p->shape() -1;
376 0 : for (Int j=0; j < nImPol_p; ++j){
377 0 : blc(polAxis_p)=j;
378 0 : trc(polAxis_p)=j;
379 0 : Slicer sl(blc, trc, Slicer::endIsLast);
380 0 : Matrix<Float> dirtySub ;
381 0 : Array<Float> buf;
382 0 : dirty_p->getSlice(buf, sl, true);
383 0 : dirtySub.reference(buf);
384 0 : matClean_p.setDirty(dirtySub);
385 0 : matClean_p.makeDirtyScales();
386 0 : Array<Float> bufMod;
387 0 : getModel=modelimage.getSlice(bufMod, sl, true);
388 0 : subModel.reference(bufMod);
389 0 : if(mask_p !=0){
390 0 : if((nMaskPol_p >1)){
391 0 : Matrix<Float> maskSub;
392 0 : Array<Float> buf2;
393 0 : mask_p->getSlice(buf2, sl, true);
394 0 : maskSub.reference(buf2);
395 0 : matClean_p.setMask(maskSub, maskThresh_p);
396 : }
397 0 : matClean_p.makeScaleMasks();
398 : }
399 0 : result=matClean_p.clean(subModel, true);
400 : //Update the private flux and residuals here
401 0 : maxResidual_p=max(maxResidual_p, matClean_p.strengthOptimum());
402 0 : if(!getModel)
403 0 : modelimage.putSlice(bufMod, sl.start());
404 0 : }
405 : }
406 : //1 pol or less
407 : else{
408 : //psf and mask should have been set above if
409 0 : Matrix<Float> dirtySub ;
410 0 : Array<Float> buf;
411 0 : dirty_p->get(buf, true);
412 0 : if(buf.shape().nelements() != 2){
413 0 : throw(AipsError("Non-expected axes in this image"));
414 : }
415 0 : dirtySub.reference(buf);
416 0 : matClean_p.setDirty(dirtySub);
417 0 : matClean_p.makeDirtyScales();
418 0 : Matrix<Float> subModel;
419 0 : Array<Float> buf0;
420 0 : Bool getModel=false;
421 0 : getModel=modelimage.get(buf0, true);
422 0 : subModel.reference(buf0);
423 0 : result=matClean_p.clean(subModel, true);
424 : //Update the private flux and residuals here
425 0 : maxResidual_p=max(maxResidual_p, matClean_p.strengthOptimum());
426 0 : if(!getModel)
427 0 : modelimage.put(subModel);
428 : }
429 0 : return result;
430 :
431 : }
432 :
433 :
434 : } //# NAMESPACE CASA - END
|
