Line data Source code
1 : //# CleanImageSkyModel.cc: Implementation of CleanImageSkyModel classes
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <casacore/casa/Arrays/ArrayMath.h>
29 : #include <synthesis/MeasurementComponents/CleanImageSkyModel.h>
30 : #include <components/ComponentModels/SkyComponent.h>
31 : #include <components/ComponentModels/ComponentList.h>
32 : #include <casacore/casa/OS/File.h>
33 : #include <synthesis/MeasurementEquations/SkyEquation.h>
34 : #include <synthesis/TransformMachines/StokesImageUtil.h>
35 : #include <casacore/casa/Exceptions/Error.h>
36 : #include <casacore/casa/BasicSL/String.h>
37 : #include <casacore/casa/Utilities/Assert.h>
38 :
39 : #include <casacore/tables/Tables/TableLock.h>
40 :
41 : #include <sstream>
42 :
43 : #include <casacore/casa/Logging/LogMessage.h>
44 : #include <casacore/casa/Logging/LogIO.h>
45 : #include <casacore/casa/Logging/LogSink.h>
46 :
47 : using namespace casacore;
48 : namespace casa { //# NAMESPACE CASA - BEGIN
49 : #define NEED_UNDERSCORES
50 : #if defined(NEED_UNDERSCORES)
51 : #define maximg maximg_
52 : #endif
53 : extern "C" {
54 : void maximg(Float*, int*, Float*, int*, int*, int*, Float*, Float*);
55 : };
56 :
57 0 : CleanImageSkyModel::CleanImageSkyModel() : ImageSkyModel(), doPolJoint_p(true)
58 : {
59 :
60 :
61 0 : }
62 :
63 0 : Bool CleanImageSkyModel::addMask(Int thismodel, ImageInterface<Float>& mask)
64 : {
65 0 : LogIO os(LogOrigin("CleanImageSkyModel", "addMask"));
66 0 : if(thismodel>=nmodels_p||thismodel<0) {
67 0 : os << LogIO::SEVERE << "Illegal model slot" << thismodel << LogIO::POST;
68 0 : return false;
69 : }
70 0 : if(Int(mask_p.nelements())<=thismodel) mask_p.resize(thismodel+1);
71 0 : mask_p[thismodel] = &mask;
72 0 : AlwaysAssert(mask_p[thismodel], AipsError);
73 0 : return true;
74 0 : }
75 :
76 0 : CleanImageSkyModel::CleanImageSkyModel(const CleanImageSkyModel& other) : ImageSkyModel() {
77 0 : operator=(other);
78 0 : };
79 :
80 0 : CleanImageSkyModel::~CleanImageSkyModel() {
81 0 : };
82 :
83 0 : CleanImageSkyModel& CleanImageSkyModel::operator=(const CleanImageSkyModel& other) {
84 0 : if(this!=&other) {
85 0 : for (Int thismodel=0;thismodel<nmodels_p;thismodel++) {
86 0 : mask_p[thismodel]=other.mask_p[thismodel];
87 0 : fluxmask_p[thismodel]=other.fluxmask_p[thismodel];
88 : }
89 : };
90 0 : return *this;
91 : }
92 :
93 0 : Bool CleanImageSkyModel::add(ComponentList& compList)
94 : {
95 0 : return ImageSkyModel::add(compList);
96 : }
97 :
98 0 : Int CleanImageSkyModel::add(ImageInterface<Float>& image, const Int maxNumXfr)
99 : {
100 0 : Int index=ImageSkyModel::add(image, maxNumXfr);
101 :
102 0 : mask_p.resize(nmodels_p);
103 0 : fluxmask_p.resize(nmodels_p);
104 :
105 0 : mask_p[index]=0;
106 0 : fluxmask_p[index]=0;
107 :
108 0 : return index;
109 : }
110 :
111 0 : Bool CleanImageSkyModel::hasMask(Int model)
112 : {
113 0 : if(mask_p.nelements()==0) return false;
114 0 : return (mask_p[model]);
115 : }
116 :
117 0 : ImageInterface<Float>& CleanImageSkyModel::mask(Int model)
118 : {
119 0 : AlwaysAssert(nmodels_p>0, AipsError);
120 0 : AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
121 0 : AlwaysAssert(mask_p[model], AipsError);
122 0 : return *mask_p[model];
123 : };
124 :
125 0 : Bool CleanImageSkyModel::hasFluxMask(Int model)
126 : {
127 0 : if(fluxmask_p.nelements()==0) return false;
128 0 : return (fluxmask_p[model]);
129 : }
130 :
131 0 : ImageInterface<Float>& CleanImageSkyModel::fluxMask(Int model)
132 : {
133 0 : AlwaysAssert(nmodels_p>0, AipsError);
134 0 : AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
135 0 : AlwaysAssert(fluxmask_p[model], AipsError);
136 0 : return *fluxmask_p[model];
137 : };
138 :
139 0 : Bool CleanImageSkyModel::addFluxMask(Int thismodel, ImageInterface<Float>& fluxMask)
140 : {
141 0 : LogIO os(LogOrigin("CleanImageSkyModel", "add"));
142 0 : if(thismodel>=nmodels_p||thismodel<0) {
143 0 : os << LogIO::SEVERE << "Illegal model slot" << thismodel << LogIO::POST;
144 0 : return false;
145 : }
146 0 : if(Int(fluxmask_p.nelements())<=thismodel) fluxmask_p.resize(thismodel+1);
147 0 : fluxmask_p[thismodel] = &fluxMask;
148 0 : AlwaysAssert(fluxmask_p[thismodel], AipsError);
149 0 : return true;
150 0 : }
151 :
152 0 : void CleanImageSkyModel::setJointStokesClean(Bool joint) {
153 0 : doPolJoint_p=joint;
154 0 : }
155 :
156 : // Find maximum residual
157 0 : Float CleanImageSkyModel::maxField(Vector<Float>& imagemax,
158 : Vector<Float>& imagemin) {
159 :
160 0 : LogIO os(LogOrigin("ImageSkyModel","maxField"));
161 :
162 0 : Float absmax=0.0;
163 0 : imagemax.resize(numberOfModels());
164 0 : imagemin.resize(numberOfModels());
165 0 : imagemax=-1e20;
166 0 : imagemin=1e20;
167 :
168 :
169 :
170 :
171 : // Loop over all models
172 0 : for (Int model=0;model<numberOfModels();model++) {
173 : // Find maximum of ggS for scaling in maximg
174 0 : Float maxggS=0.0;
175 : {
176 0 : LatticeExprNode LEN = max(LatticeExpr<Float>(ggS(model)));
177 0 : maxggS = LEN.getFloat();
178 0 : }
179 : // Remember that the residual image can be either as specified
180 : // or created specially.
181 0 : ImageInterface<Float>* imagePtr=0;
182 0 : if(residual_p[model]) {
183 0 : imagePtr=residual_p[model];
184 : }
185 : else {
186 0 : imagePtr=(ImageInterface<Float> *)residualImage_p[model];
187 : }
188 0 : AlwaysAssert(imagePtr, AipsError);
189 0 : AlwaysAssert(imagePtr->shape().nelements()==4, AipsError);
190 0 : Int nx=imagePtr->shape()(0);
191 0 : Int ny=imagePtr->shape()(1);
192 0 : Int npol=imagePtr->shape()(2);
193 0 : Int nchan=imagePtr->shape()(3);
194 :
195 : // AlwaysAssert((npol==1)||(npol==2)||(npol==4), AipsError);
196 :
197 : // Loop over all channels
198 0 : IPosition oneSlab(4, nx, ny, npol, 1);
199 0 : LatticeStepper ls(imagePtr->shape(), oneSlab, IPosition(4, 0, 1, 2, 3));
200 0 : IPosition onePlane(4, nx, ny, 1, 1);
201 0 : LatticeStepper fsls(ggS(model).shape(), oneSlab,
202 0 : IPosition(4,0,1,2,3));
203 0 : RO_LatticeIterator<Float> ggSli(ggS(model),fsls);
204 0 : LatticeIterator<Float> imageli(*imagePtr, ls);
205 :
206 : // If we are using a mask then reset the region to be
207 : // cleaned
208 0 : Array<Float> maskArray;
209 0 : RO_LatticeIterator<Float> maskIter;
210 0 : Bool cubeMask=false;
211 :
212 0 : Int domask=0;
213 0 : if(hasMask(model)) {
214 0 : Int mx=mask(model).shape()(0);
215 0 : Int my=mask(model).shape()(1);
216 0 : Int mpol=mask(model).shape()(2);
217 0 : Int nMaskChan=0;
218 0 : if(mask(model).shape().nelements()==4){
219 0 : nMaskChan=mask(model).shape()(3);
220 : }
221 0 : if( (nchan >1) && nMaskChan==nchan)
222 0 : cubeMask=true;
223 0 : if((mx != nx) || (my != ny) || (mpol != npol)){
224 0 : throw(AipsError("Mask image shape is not the same as dirty image"));
225 : }
226 0 : LatticeStepper mls(mask(model).shape(), oneSlab,
227 0 : IPosition(4, 0, 1, 2, 3));
228 :
229 0 : RO_LatticeIterator<Float> maskli(mask(model), mls);
230 0 : maskli.reset();
231 0 : maskIter=maskli;
232 0 : if (maskli.cursor().shape().nelements() > 1) {
233 0 : domask=1;
234 0 : maskArray=maskli.cursor();
235 : }
236 0 : }
237 :
238 0 : Int chan=0;
239 : Float imax, imin;
240 0 : imax=-1E20; imagemax(model)=imax;
241 0 : imin=+1E20; imagemin(model)=imin;
242 :
243 0 : for (imageli.reset(),ggSli.reset();!imageli.atEnd();imageli++,ggSli++,chan++) {
244 : Float fmax, fmin;
245 : Bool delete_its;
246 : Bool delete_its2;
247 : // Renormalize by the weights
248 0 : if(cubeMask){
249 0 : if(maskIter.cursor().shape().nelements() > 1){
250 0 : domask=1;
251 0 : maskArray=maskIter.cursor();
252 : }
253 0 : maskIter++;
254 :
255 : }
256 :
257 0 : Cube<Float> weight(sqrt(ggSli.cursor().nonDegenerate(3)/maxggS));
258 : // resid=(imageli.cursor().nonDegenerate(3)) does not make a
259 : // copy. This is a bug in LatticeIterator.cursor().
260 : //
261 : // As a result resid gets modified for the rest of the code
262 : // (residual image is multipled by sqrt(ggS)!!
263 : //
264 : // For now, using Cube::assign() to force a copy.
265 0 : Cube<Float>resid;resid.assign(imageli.cursor().nonDegenerate(3));
266 0 : for (Int pol=0;pol<npol;pol++) {
267 0 : for (Int iy=0;iy<ny;iy++) {
268 0 : for (Int ix=0;ix<nx;ix++) {
269 0 : resid(ix,iy,pol)*=weight(ix,iy,pol);
270 : }
271 : }
272 : }
273 0 : const Float* limage_data=resid.getStorage(delete_its);
274 0 : const Float* lmask_data=maskArray.getStorage(delete_its2);
275 0 : maximg((Float*)limage_data, &domask, (Float*)lmask_data,
276 : &nx, &ny, &npol, &fmin, &fmax);
277 :
278 :
279 0 : resid.freeStorage(limage_data, delete_its);
280 0 : maskArray.freeStorage(lmask_data, delete_its2);
281 0 : if(fmax<0.99*imax) fmax=0.0;
282 0 : if(fmin>0.99*imin) fmin=0.0;
283 0 : if(abs(fmax)>absmax) absmax=abs(fmax);
284 0 : if(abs(fmin)>absmax) absmax=abs(fmin);
285 0 : if(fmin<imagemin(model)) imagemin(model)=fmin;
286 0 : if(fmax>imagemax(model)) imagemax(model)=fmax;
287 0 : }
288 0 : }
289 0 : return absmax;
290 0 : };
291 :
292 :
293 : } //# NAMESPACE CASA - END
294 :
|
