Line data Source code
1 : //# MSCleanImageSkyModel.cc: Implementation of MSCleanImageSkyModel class
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,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: aips2-request@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 :
29 : #include <casacore/casa/Arrays/ArrayMath.h>
30 : #include <synthesis/MeasurementComponents/MSCleanImageSkyModel.h>
31 : #include <casacore/casa/OS/File.h>
32 : #include <synthesis/MeasurementEquations/ImageMSCleaner.h>
33 : #include <casacore/images/Images/SubImage.h>
34 : #include <casacore/lattices/LRegions/LCBox.h>
35 : #include <synthesis/MeasurementEquations/SkyEquation.h>
36 : #include <synthesis/MeasurementEquations/LatticeModel.h>
37 : #include <synthesis/MeasurementEquations/LatConvEquation.h>
38 : #include <casacore/lattices/LEL/LatticeExprNode.h>
39 : #include <casacore/casa/Exceptions/Error.h>
40 : #include <casacore/casa/BasicSL/String.h>
41 : #include <casacore/casa/Utilities/Assert.h>
42 : #include <casacore/lattices/Lattices/TempLattice.h>
43 :
44 : #include <sstream>
45 :
46 : #include <casacore/casa/Logging/LogMessage.h>
47 : #include <casacore/casa/Logging/LogSink.h>
48 :
49 :
50 : using namespace casacore;
51 : namespace casa { //# NAMESPACE CASA - BEGIN
52 :
53 0 : MSCleanImageSkyModel::MSCleanImageSkyModel(const Int nscales, const Int stoplargenegatives, const Int stoppointmode, const Float smallScaleBias)
54 0 : : method_p(NSCALES), nscales_p(nscales), userScaleSizes_p(0), stopLargeNegatives_p(stoplargenegatives), stopPointMode_p(stoppointmode), smallScaleBias_p(smallScaleBias)
55 : {
56 0 : modified_p=true;
57 0 : donePSF_p=false;
58 :
59 0 : };
60 :
61 0 : MSCleanImageSkyModel::MSCleanImageSkyModel(const Vector<Float>& userScaleSizes, const Int stoplarge, const Int stoppoint, const Float smallScaleBias)
62 : : method_p(USERVECTOR), nscales_p(0), userScaleSizes_p(userScaleSizes),
63 : stopLargeNegatives_p(stoplarge), stopPointMode_p(stoppoint),
64 0 : smallScaleBias_p(smallScaleBias)
65 : {
66 0 : modified_p=true;
67 0 : donePSF_p=false;
68 :
69 0 : };
70 :
71 0 : MSCleanImageSkyModel::~MSCleanImageSkyModel()
72 : {
73 :
74 0 : };
75 :
76 : // Clean solver
77 0 : Bool MSCleanImageSkyModel::solve(SkyEquation& se) {
78 :
79 :
80 0 : LogIO os(LogOrigin("MSCleanImageSkyModel","solve"));
81 :
82 0 : if(numberOfModels()>1) {
83 0 : os << "Cannot process more than one field" << LogIO::EXCEPTION;
84 : }
85 : /*
86 : if (displayProgress_p) {
87 : progress_p = new LatticeCleanProgress( pgplotter_p );
88 : }
89 : */
90 : // Make the residual image
91 0 : if(modified_p)
92 0 : makeNewtonRaphsonStep(se, false, (numberIterations()<1)?true:False);
93 :
94 : //Make the PSF
95 0 : if(!donePSF_p)
96 0 : makeApproxPSFs(se);
97 :
98 :
99 0 : if(numberIterations() <1){
100 0 : return true;
101 : }
102 :
103 :
104 0 : if(!isSolveable(0)) {
105 0 : os << "Model 1 is not solveable!" << LogIO::EXCEPTION;
106 : }
107 :
108 :
109 0 : Vector<Float> scaleSizes(0);
110 0 : if (method_p == USERVECTOR) {
111 0 : if (userScaleSizes_p.nelements() <= 0) {
112 : os << LogIO::SEVERE
113 : << "Need at least one scale for method uservector"
114 0 : << LogIO::POST;
115 : }
116 0 : os << "Creating scales from uservector method: " << LogIO::POST;
117 0 : for(uInt scale=0; scale < userScaleSizes_p.nelements(); scale++) {
118 0 : os << "scale " << scale << " = " << userScaleSizes_p(scale)
119 0 : << " pixels" << LogIO::POST;
120 : }
121 : } else {
122 0 : if (nscales_p <= 0) nscales_p = 1;
123 0 : scaleSizes.resize(nscales_p);
124 : os << "Creating " << nscales_p <<
125 0 : " scales from powerlaw nscales method" << LogIO::POST;
126 0 : scaleSizes(0) = 0.0;
127 0 : os << "scale 1 = 0.0 pixels " << LogIO::POST;
128 0 : Float scaleInc = 2.0;
129 0 : for (Int scale=1; scale<nscales_p;scale++) {
130 0 : scaleSizes(scale) =
131 0 : scaleInc * pow(10.0, (Float(scale)-2.0)/2.0);
132 0 : os << "scale " << scale+1 << " = " << scaleSizes(scale)
133 0 : << " pixels" << LogIO::POST;
134 : }
135 : }
136 :
137 0 : Int npol=image(0).shape()(2);
138 0 : Int nchan=image(0).shape()(3);
139 :
140 0 : AlwaysAssert((npol==1)||(npol==2)||(npol==3)||(npol==4), AipsError);
141 0 : IPosition blcDirty(image(0).shape().nelements(), 0);
142 0 : IPosition trcDirty(image(0).shape()-1);
143 :
144 0 : if(hasMask(0)) {
145 0 : Int masknpol=mask(0).shape()(2);
146 0 : if(masknpol>1) {
147 0 : if(masknpol!=npol) {
148 0 : os << "Mask has more than one polarization but not the same as the image" << LogIO::EXCEPTION;
149 : }
150 : else {
151 0 : os << "Mask is a cube in polarization - will use appropriate plane for each polarization" << LogIO::POST;
152 : }
153 : }
154 0 : Int masknchan=mask(0).shape()(3);
155 0 : if(masknchan>1) {
156 0 : if(masknchan!=nchan) {
157 0 : os << "Mask has more than one channel but not the same as the image" << LogIO::EXCEPTION;
158 : }
159 : else {
160 0 : os << "Mask is a spectral cube - will use appropriate plane for each channel" << LogIO::POST;
161 : }
162 : }
163 : }
164 :
165 0 : Int iterUsed=0;
166 0 : Float maxRes=0.0;
167 0 : Int converged=0;
168 : // Loop over all channels and polarizations
169 0 : for (Int chan=0; chan<nchan; chan++) {
170 :
171 0 : if(nchan>1) {
172 0 : os<<"Processing channel "<<chan<<" of 0 to "<<nchan-1<<LogIO::POST;
173 : }
174 :
175 : // Load the PSF for this channel
176 : // Decide whether to clean this channel or not, based on PSF peak being non-zero.
177 0 : blcDirty(3) = chan; trcDirty(3) = chan;
178 0 : blcDirty(2) = 0; trcDirty(2) = 0;
179 0 : LCBox firstPolPlane(blcDirty, trcDirty, image(0).shape());
180 0 : SubImage<Float> subPsf(PSF(0), firstPolPlane);
181 : Float psfmax;
182 : {
183 0 : LatticeExprNode node = max(subPsf);
184 0 : psfmax = node.getFloat();
185 : }
186 :
187 0 : if(psfmax==0.0) { // PSF is not valid
188 0 : os << "No data for this channel: skipping" << LogIO::POST;
189 : } else { // PSF is valid.
190 0 : SubImage<Float> subMask;
191 0 : for (Int pol=0; pol<npol; pol++) { // Run MS-Clean on each polarization
192 0 : blcDirty(2) = pol; trcDirty(2) = pol;
193 0 : if(npol>1) {
194 0 : os<<"Processing polarization "<<pol+1<<" of "<<npol<<LogIO::POST;
195 : }
196 :
197 :
198 0 : LCBox onePlane(blcDirty, trcDirty, image(0).shape());
199 0 : SubImage<Float> subDirty(residual(0), onePlane, true);
200 0 : ImageMSCleaner cleaner(subPsf, subDirty);
201 :
202 :
203 0 : Bool doClean=true;
204 0 : String algorithm="msclean";
205 0 : if(hasMask(0)) {
206 0 : IPosition blcMask(mask(0).shape().nelements(), 0);
207 0 : IPosition trcMask(mask(0).shape()-1);
208 :
209 0 : Int masknpol=mask(0).shape()(2);
210 0 : Int masknchan=mask(0).shape()(3);
211 :
212 0 : if(masknpol==npol) {
213 0 : blcMask(2)=pol;
214 0 : trcMask(2)=pol;
215 : }
216 : else {
217 0 : blcMask(2)=0;
218 0 : trcMask(2)=0;
219 : }
220 0 : if(masknchan==nchan) {
221 0 : blcMask(3)=chan;
222 0 : trcMask(3)=chan;
223 : }
224 : else {
225 0 : blcMask(3)=0;
226 0 : trcMask(3)=0;
227 : }
228 0 : LCBox maskPlane(blcMask, trcMask, mask(0).shape());
229 0 : subMask=SubImage<Float>( mask(0), maskPlane, false);
230 : // Check for empty mask
231 0 : LatticeExprNode sumMask = sum(subMask);
232 0 : if(sumMask.getFloat()==0.0) {
233 0 : os << LogIO::WARN << "Mask is specified but empty - no cleaning performed" << LogIO::POST;
234 0 : doClean=false;
235 : }
236 : else {
237 0 : cleaner.setMask(subMask);
238 : //Using mask so the user knows best.
239 0 : cleaner.ignoreCenterBox(true);
240 0 : algorithm="fullmsclean";
241 : }
242 : }
243 :
244 :
245 :
246 0 : if(doClean) {
247 0 : SubImage<Float> subImage(image(0), onePlane, true);
248 :
249 0 : if (method_p == USERVECTOR) {
250 0 : cleaner.setscales(userScaleSizes_p);
251 : } else {
252 0 : cleaner.setscales(scaleSizes);
253 : }
254 0 : cleaner.setSmallScaleBias(smallScaleBias_p);
255 0 : cleaner.stopPointMode(stopPointMode_p);
256 : //cleaner.setcontrol(CleanEnums::MULTISCALE, numberIterations(), gain(),
257 : // Quantity(threshold(), "Jy"), true);
258 0 : if(stopLargeNegatives_p >0)
259 0 : cleaner.stopAtLargeScaleNegative();
260 :
261 0 : converged=cleaner.clean(subImage, algorithm, numberIterations(),gain(), Quantity(threshold(), "Jy"), Quantity(0.0, "%"), true);
262 0 : Int stoplarge=stopLargeNegatives_p;
263 0 : while( (converged==-2) && stoplarge > 0){
264 0 : --stoplarge;
265 0 : converged=cleaner.clean(subImage, algorithm, numberIterations(),gain(),
266 0 : Quantity(threshold(), "Jy"), Quantity(0, "%"),true);
267 : }
268 : // calculate residuals
269 :
270 0 : LatticeModel lm(subImage);
271 0 : LatConvEquation eqn(subPsf, subDirty);
272 0 : TempLattice<Float> restl( subDirty.shape() );
273 0 : eqn.residual(restl, lm);
274 0 : subDirty.copyData(restl);
275 0 : iterUsed=max(iterUsed, cleaner.iteration());
276 : }
277 : }// end of polarization loop
278 : }// end of if (valid psf)
279 : }// end of channel loop
280 0 : Vector<Float> minres;
281 0 : Vector <Float> maxres;
282 0 : maxRes=maxField(maxres, minres);
283 0 : setThreshold(maxRes);
284 :
285 0 : modified_p=true;
286 :
287 :
288 0 : return(converged);
289 : };
290 :
291 : } //# NAMESPACE CASA - END
292 :
|
