Line data Source code
1 : //# Copyright (C) 2015
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 General Public License as published by
6 : //# the Free Software Foundation; either version 3 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/MatrixNACleaner.h>
27 : #include <synthesis/MeasurementEquations/ImageNACleaner.h>
28 : #include <casacore/images/Images/ImageInterface.h>
29 : #include <casacore/images/Images/PagedImage.h>
30 : #include <casacore/images/Images/TempImage.h>
31 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
32 : #include <casacore/coordinates/Coordinates/CoordinateUtil.h>
33 : using namespace std;
34 : using namespace casacore;
35 : namespace casa { //# NAMESPACE CASA - BEGIN
36 :
37 0 : ImageNACleaner::ImageNACleaner(): psf_p(nullptr), dirty_p(nullptr), mask_p(nullptr), nPsfChan_p(0),
38 : nImChan_p(0), nPsfPol_p(0), nImPol_p(0), chanAxis_p(-1),
39 0 : polAxis_p(-1), nMaskChan_p(0), nMaskPol_p(0), maxResidual_p(0.0)
40 : {
41 :
42 :
43 0 : }
44 0 : ImageNACleaner::ImageNACleaner(ImageInterface<Float>& psf,
45 0 : ImageInterface<Float>& dirty): mask_p(nullptr),
46 0 : nMaskChan_p(0), nMaskPol_p(0),maxResidual_p(0.0){
47 0 : psf_p=CountedPtr<ImageInterface<Float> >(&psf, false);
48 0 : dirty_p=CountedPtr<ImageInterface<Float> >(&dirty, false);
49 0 : chanAxis_p=CoordinateUtil::findSpectralAxis(dirty_p->coordinates());
50 0 : Vector<Stokes::StokesTypes> whichPols;
51 0 : polAxis_p=CoordinateUtil::findStokesAxis(whichPols, dirty_p->coordinates());
52 0 : if(chanAxis_p > -1){
53 0 : nPsfChan_p=psf_p->shape()(chanAxis_p);
54 0 : nImChan_p=dirty_p->shape()(chanAxis_p);
55 : }
56 0 : if(polAxis_p > -1){
57 0 : nPsfPol_p=psf_p->shape()(polAxis_p);
58 0 : nImPol_p=dirty_p->shape()(polAxis_p);
59 : }
60 0 : }
61 :
62 :
63 0 : ImageNACleaner::ImageNACleaner(const ImageNACleaner& other){
64 0 : operator=(other);
65 0 : }
66 :
67 0 : ImageNACleaner::~ImageNACleaner(){
68 :
69 0 : }
70 :
71 0 : ImageNACleaner& ImageNACleaner::operator=(const ImageNACleaner& other){
72 0 : if (this != &other) {
73 0 : matClean_p=other.matClean_p;
74 0 : psf_p=other.psf_p;
75 0 : dirty_p=other.dirty_p;
76 0 : mask_p=other.mask_p;
77 0 : nPsfChan_p=other.nPsfChan_p;
78 0 : nImChan_p=other.nImChan_p;
79 0 : nPsfPol_p=other.nPsfPol_p;
80 0 : nImPol_p=other.nImPol_p;
81 0 : chanAxis_p=other.chanAxis_p;
82 0 : nMaskChan_p=other.nMaskChan_p;
83 0 : nMaskPol_p=other.nMaskPol_p;
84 0 : polAxis_p=other.polAxis_p;
85 0 : maxResidual_p=other.maxResidual_p;
86 : }
87 0 : return *this;
88 : }
89 0 : void ImageNACleaner::setDirty(ImageInterface<Float>& dirty){
90 0 : dirty_p=CountedPtr<ImageInterface<Float> >(&dirty, false);
91 0 : chanAxis_p=CoordinateUtil::findSpectralAxis(dirty_p->coordinates());
92 0 : Vector<Stokes::StokesTypes> whichPols;
93 0 : polAxis_p=CoordinateUtil::findStokesAxis(whichPols, dirty_p->coordinates());
94 0 : if(chanAxis_p > -1){
95 0 : nImChan_p=dirty_p->shape()(chanAxis_p);
96 : }
97 : else
98 0 : nImChan_p=0;
99 0 : if(polAxis_p > -1){
100 0 : nImPol_p=dirty_p->shape()(polAxis_p);
101 : }
102 : else
103 0 : nImPol_p=0;
104 0 : }
105 :
106 0 : void ImageNACleaner::setPsf(ImageInterface<Float> & psf){
107 0 : psf_p=CountedPtr<ImageInterface<Float> >(&psf,false);
108 0 : Int chanAxis=CoordinateUtil::findSpectralAxis(psf_p->coordinates());
109 0 : Vector<Stokes::StokesTypes> whichPols;
110 0 : Int polAxis=CoordinateUtil::findStokesAxis(whichPols, psf_p->coordinates());
111 0 : if(chanAxis > -1){
112 0 : nPsfChan_p=psf_p->shape()(chanAxis);
113 : }
114 : else
115 0 : nPsfChan_p=0;
116 0 : if(polAxis > -1){
117 0 : nPsfPol_p=psf_p->shape()(polAxis);
118 : }
119 : else
120 0 : nPsfPol_p=0;
121 0 : }
122 :
123 :
124 0 : void ImageNACleaner::setcontrol(const Int niter,
125 : const Float gain, const Quantity& aThreshold,
126 : const Int supp, const Int memtype, const Float numsigma){
127 0 : matClean_p.setcontrol(niter, gain, aThreshold, supp, memtype, numsigma);
128 0 : }
129 :
130 0 : Int ImageNACleaner::iteration() const{
131 0 : return matClean_p.iteration();
132 : }
133 :
134 : /* void ImageMSCleaner::startingIteration(const Int starting){
135 : matClean_p.startingIteration(starting);
136 : }
137 : */
138 0 : void ImageNACleaner::setMask(ImageInterface<Float> & mask){
139 :
140 :
141 0 : mask_p=CountedPtr<ImageInterface<Float> >(&mask, false);
142 0 : Int chanAxis=CoordinateUtil::findSpectralAxis(mask_p->coordinates());
143 0 : Vector<Stokes::StokesTypes> whichPols;
144 0 : Int polAxis=CoordinateUtil::findStokesAxis(whichPols, mask_p->coordinates());
145 0 : if(chanAxis > -1){
146 0 : nMaskChan_p=mask_p->shape()(chanAxis);
147 : }
148 : else
149 0 : nMaskChan_p=0;
150 0 : if(polAxis > -1){
151 0 : nMaskPol_p=mask_p->shape()(polAxis);
152 : }
153 : else
154 0 : nMaskPol_p=0;
155 :
156 :
157 0 : }
158 :
159 :
160 :
161 0 : Bool ImageNACleaner::setupMatCleaner(const Int niter,
162 : const Float gain, const Quantity& threshold, const Int masksupp, const Int memType, const Float numsigma){
163 :
164 0 : LogIO os(LogOrigin("ImageNACleaner", "setupMatCleaner()", WHERE));
165 :
166 0 : matClean_p.setcontrol(niter, gain, threshold, masksupp, memType, numsigma);
167 :
168 :
169 0 : return true;
170 : }
171 :
172 0 : Int ImageNACleaner::clean(ImageInterface<Float> & modelimage,
173 : const Int niter,
174 : const Float gain, const Quantity& threshold, const Int masksupp, const Int memType, const Float numSigma){
175 :
176 :
177 :
178 :
179 :
180 0 : Int result=0;
181 : ///case of single plane mask
182 : //now may be a time to set stuff scales will be done later
183 0 : if(!setupMatCleaner(niter, gain, threshold, masksupp, memType, numSigma))
184 0 : return false;
185 : //cerr << "nPol " << nMaskPol_p << " " << nPsfPol_p << " " << nImPol_p << endl;
186 : //cerr << "nChan " << nMaskChan_p << " " << nPsfChan_p << " " << nImChan_p << endl;
187 0 : if( (nMaskPol_p >1) && (nMaskPol_p != nImPol_p))
188 0 : throw(AipsError("Donot know how to deal with mask that has different pol axis as Image"));
189 0 : CountedPtr<ImageInterface<Float> > somemask;
190 0 : if(!mask_p){
191 : //make one
192 0 : somemask=CountedPtr<TempImage<Float> > (new TempImage<Float>((*dirty_p).shape(), (*dirty_p).coordinates())) ;
193 0 : setMask(*somemask);
194 : }
195 :
196 : /////TEST
197 : //PagedImage<Float> someResid((*dirty_p).shape(), (*dirty_p).coordinates(), "decon.resid");
198 : /////
199 0 : if(!psf_p)
200 0 : throw(AipsError("No PSF defined "));
201 : ///case of single plane psf
202 0 : if(nPsfChan_p < 2){
203 0 : Matrix<Float> psfMat;
204 0 : Array<Float> mbuf;
205 0 : if(nPsfPol_p > 1){
206 : ///First stokes psf is good enough
207 0 : IPosition blc(dirty_p->ndim(),0);
208 0 : IPosition trc=(dirty_p->shape()) -1;
209 0 : trc(polAxis_p)=0;
210 0 : Slicer sl(blc, trc, Slicer::endIsLast);
211 0 : psf_p->getSlice(mbuf, sl, true);
212 : }
213 : else{
214 0 : psf_p->get(mbuf, true);
215 : }
216 0 : psfMat.reference(mbuf);
217 0 : matClean_p.setPsf(psfMat);
218 :
219 : }
220 :
221 :
222 :
223 : //has cube axis
224 0 : if(chanAxis_p > -1 && nImChan_p > 0){
225 0 : Int npol=1;
226 0 : if(polAxis_p > -1 && nImPol_p >0){
227 0 : npol=nImPol_p;
228 : }
229 0 : Matrix<Float> subModel;
230 0 : Bool getModel=false;
231 0 : Bool getMask=false;
232 0 : IPosition blc(dirty_p->ndim(), 0);
233 0 : IPosition trc=dirty_p->shape() -1;
234 0 : for (Int k=0; k < nImChan_p ; ++k){
235 0 : for (Int j=0; j < npol; ++j){
236 0 : setupMatCleaner(niter, gain, threshold, masksupp, memType, numSigma);
237 0 : if(npol > 1 ){
238 0 : blc(polAxis_p)=j;
239 0 : trc(polAxis_p)=j;
240 : }
241 0 : blc(chanAxis_p)=k;
242 0 : trc(chanAxis_p)=k;
243 : //cerr << "blc " << blc << " trc " << trc << endl;
244 0 : Slicer sl(blc, trc, Slicer::endIsLast);
245 0 : Matrix<Float> dirtySub ;
246 0 : Array<Float> buf;
247 0 : dirty_p->getSlice(buf, sl, true);
248 0 : if(dirty_p->isMasked()){
249 0 : cerr << "zeroing masked dirty" << endl;
250 0 : Array<Bool> bitmask;
251 0 : dirty_p->pixelMask().getSlice(bitmask, sl, true);
252 0 : buf(!bitmask)=0.0;
253 0 : matClean_p.setPixFlag(bitmask);
254 : }
255 0 : dirtySub.reference(buf);
256 0 : matClean_p.setDirty(dirtySub);
257 0 : Array<Float> bufMod;
258 0 : getModel=modelimage.getSlice(bufMod, sl, true);
259 0 : subModel.reference(bufMod);
260 :
261 0 : if((nPsfChan_p>1) && (k<(nPsfChan_p-1))){
262 0 : Matrix<Float> psfSub;
263 0 : Array<Float> buf1;
264 0 : psf_p->getSlice(buf1, sl, true);
265 0 : psfSub.reference(buf1);
266 0 : matClean_p.setPsf(psfSub);
267 : }
268 0 : Array<Float> buf2;
269 0 : if(mask_p){
270 : // if((nMaskChan_p >1) && (k<(nMaskChan_p))){
271 0 : Matrix<Float> maskSub;
272 0 : getMask=mask_p->getSlice(buf2, sl, true);
273 : //cerr << "getmask " << getMask << " buf2 " << buf2.shape() << endl;
274 0 : maskSub.reference(buf2);
275 0 : matClean_p.setMask(maskSub);
276 : //Set mask above does the makeScaleMasks as psf-scale-Xfr are valid
277 : // }
278 :
279 :
280 : // else{
281 : //matClean_p.unsetMask();
282 : //}
283 : }
284 0 : result=matClean_p.clean(subModel);
285 : //Update the private flux and residuals here
286 0 : cerr << "maxResidual " << matClean_p.maxResidual() << " iteration " << matClean_p.iteration() << endl;
287 0 : maxResidual_p=max(maxResidual_p, matClean_p.maxResidual());
288 0 : if(!getModel)
289 0 : modelimage.putSlice(bufMod, sl.start());
290 : /////TESTING
291 :
292 : //someResid.putSlice(matClean_p.getResidual().reform(bufMod.shape()), sl.start());
293 :
294 : ///////
295 : // cerr << "mask " << max(buf2) << " buf2 " << buf2.shape() << " start " << sl.start() << endl;
296 0 : if(!getMask)
297 0 : mask_p->putSlice(buf2, sl.start());
298 : }
299 0 : }
300 : }
301 : //has no channel but has pol
302 0 : else if(polAxis_p > -1 && nImPol_p > 1){
303 : //The psf should have been set above before the if loop
304 : //for this case
305 0 : Matrix<Float> subModel;
306 0 : Bool getModel=false;
307 0 : Bool getMask=false;
308 0 : Array<Float> buf2;
309 0 : IPosition blc(dirty_p->ndim(), 0);
310 0 : IPosition trc=dirty_p->shape() -1;
311 0 : for (Int j=0; j < nImPol_p; ++j){
312 0 : blc(polAxis_p)=j;
313 0 : trc(polAxis_p)=j;
314 0 : Slicer sl(blc, trc, Slicer::endIsLast);
315 0 : Matrix<Float> dirtySub ;
316 0 : Array<Float> buf;
317 0 : dirty_p->getSlice(buf, sl, true);
318 0 : dirtySub.reference(buf);
319 0 : matClean_p.setDirty(dirtySub);
320 0 : Array<Float> bufMod;
321 0 : getModel=modelimage.getSlice(bufMod, sl, true);
322 0 : subModel.reference(bufMod);
323 0 : if(!mask_p){
324 0 : if((nMaskPol_p >1)){
325 0 : Matrix<Float> maskSub;
326 :
327 0 : getMask=mask_p->getSlice(buf2, sl, true);
328 0 : maskSub.reference(buf2);
329 0 : matClean_p.setMask(maskSub);
330 : }
331 :
332 : }
333 0 : result=matClean_p.clean(subModel);
334 : //Update the private flux and residuals here
335 0 : maxResidual_p=max(maxResidual_p, matClean_p.maxResidual());
336 0 : if(!getModel)
337 0 : modelimage.putSlice(bufMod, sl.start());
338 0 : if(!getMask)
339 0 : mask_p->putSlice(buf2, sl.start());
340 0 : }
341 : }
342 : //1 pol or less
343 : else{
344 : ////Will look at this later
345 : //psf and mask should have been set above if
346 0 : Matrix<Float> dirtySub ;
347 0 : Array<Float> buf;
348 0 : dirty_p->get(buf, true);
349 0 : if(buf.shape().nelements() != 2){
350 0 : throw(AipsError("Non-expected axes in this image"));
351 : }
352 0 : dirtySub.reference(buf);
353 0 : matClean_p.setDirty(dirtySub);
354 0 : Matrix<Float> subModel;
355 0 : Array<Float> buf0;
356 0 : Bool getModel=false;
357 0 : getModel=modelimage.get(buf0, true);
358 0 : subModel.reference(buf0);
359 0 : result=matClean_p.clean(subModel);
360 : //Update the private flux and residuals here
361 0 : maxResidual_p=max(maxResidual_p, matClean_p.maxResidual());
362 0 : if(!getModel)
363 0 : modelimage.put(subModel);
364 : }
365 0 : return result;
366 :
367 : }
368 :
369 :
370 : } //# NAMESPACE CASA - END
|
