Line data Source code
1 : //# ClarkCleanImageSkyModel.cc: Implementation of ClarkCleanImageSkyModel 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: 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/ClarkCleanImageSkyModel.h>
30 : #include <casacore/images/Images/PagedImage.h>
31 : #include <casacore/casa/OS/File.h>
32 : #include <casacore/lattices/Lattices/LatticeStepper.h>
33 : #include <casacore/lattices/Lattices/LatticeIterator.h>
34 : #include <casacore/lattices/LEL/LatticeExpr.h>
35 : #include <casacore/lattices/LEL/LatticeExprNode.h>
36 : #include <casacore/lattices/LRegions/LCBox.h>
37 : #include <casacore/lattices/Lattices/SubLattice.h>
38 : #include <synthesis/MeasurementEquations/SkyEquation.h>
39 : #include <casacore/casa/Exceptions/Error.h>
40 : #include <casacore/casa/BasicSL/String.h>
41 : #include <casacore/casa/Utilities/Assert.h>
42 :
43 : #include <sstream>
44 :
45 : #include <casacore/casa/Logging/LogMessage.h>
46 : #include <casacore/casa/Logging/LogIO.h>
47 : #include <casacore/casa/Logging/LogSink.h>
48 :
49 : #include <synthesis/MeasurementEquations/LatConvEquation.h>
50 : #include <synthesis/MeasurementEquations/ClarkCleanLatModel.h>
51 : #include <synthesis/MeasurementEquations/ClarkCleanProgress.h>
52 :
53 :
54 : using namespace casacore;
55 : namespace casa { //# NAMESPACE CASA - BEGIN
56 :
57 0 : ClarkCleanImageSkyModel::ClarkCleanImageSkyModel() : CleanImageSkyModel(), itsProgress(0) {
58 :
59 :
60 :
61 0 : }
62 :
63 :
64 0 : ClarkCleanImageSkyModel::~ClarkCleanImageSkyModel()
65 : {
66 0 : if (itsProgress) delete itsProgress;
67 0 : };
68 :
69 : // Clean solver
70 0 : Bool ClarkCleanImageSkyModel::solve(SkyEquation& se) {
71 :
72 0 : LogIO os(LogOrigin("ClarkCleanImageSkyModel","solve",WHERE));
73 :
74 0 : Bool converged=false;
75 0 : if(numberOfModels()>1) {
76 0 : os << "Cannot process more than one field" << LogIO::EXCEPTION;
77 : }
78 :
79 : // Make the residual image
80 0 : if(modified_p)
81 0 : makeNewtonRaphsonStep(se, false, (numberIterations()<1)?true:False);
82 :
83 : //Make the PSF
84 0 : if(!donePSF_p)
85 0 : makeApproxPSFs(se);
86 :
87 0 : if(numberIterations() < 1){
88 :
89 0 : return true;
90 : }
91 :
92 0 : Float maxRes=0.0;
93 0 : Int iterused=0;
94 :
95 0 : if(hasMask(0))
96 0 : converged = clean(image(0), residual(0), PSF(0), mask(0), maxRes, iterused, gain(), numberIterations(),
97 0 : threshold(), cycleFactor_p, true, doPolJoint_p);
98 : else{
99 0 : ImageInterface<Float> *tmpMask=0;
100 0 : converged = clean(image(0), residual(0), PSF(0), *tmpMask, maxRes, iterused, gain(), numberIterations(),
101 0 : threshold(), cycleFactor_p, false, doPolJoint_p);
102 : }
103 :
104 0 : setThreshold(maxRes);
105 0 : setNumberIterations(iterused);
106 0 : modified_p=true;
107 0 : return(converged);
108 0 : };
109 :
110 :
111 0 : Bool ClarkCleanImageSkyModel::clean(ImageInterface<Float>& image, ImageInterface<Float> & residual,
112 : ImageInterface<Float>& psf, ImageInterface<Float>& mask, Float& maxresidual, Int& iterused, Float gain, Int numIter,
113 : Float thresh, Float cycleFactor, Bool useMask, Bool doPolJoint){
114 :
115 0 : Bool converged=false;
116 0 : LogIO os(LogOrigin("ClarkCleanImageSkyModel","clean",WHERE));
117 0 : Int nx=image.shape()(0);
118 0 : Int ny=image.shape()(1);
119 0 : Int npol=image.shape()(2);
120 0 : Int nchan=image.shape()(3);
121 0 : Int polloop=1;
122 0 : Vector<String> stokesID(npol, "");
123 0 : if (!doPolJoint) {
124 0 : polloop=npol;
125 0 : CoordinateSystem cs=image.coordinates();
126 0 : Int stokesindex=cs.findCoordinate(Coordinate::STOKES);
127 0 : StokesCoordinate stcoord=cs.stokesCoordinate(stokesindex);
128 0 : for (Int jj=0; jj < npol ; ++jj){
129 0 : stokesID(jj)=Stokes::name(Stokes::type(stcoord.stokes()(jj)));
130 : }
131 0 : }
132 :
133 :
134 0 : Int newNx=nx;
135 0 : Int newNy=ny;
136 :
137 : Int xbeg, xend, ybeg, yend;
138 : //default clean box
139 0 : xbeg=nx/4;
140 0 : xend=3*nx/4-1;
141 0 : ybeg=ny/4;
142 0 : yend=3*ny/4-1;
143 :
144 0 : Bool isCubeMask=false;
145 : //AlwaysAssert((npol==1)||(npol==2)||(npol==4), AipsError);
146 :
147 0 : Lattice<Float>* mask_sl = 0;
148 0 : RO_LatticeIterator<Float>* maskli = 0;
149 :
150 0 : if(useMask) {
151 : // AlwaysAssert(mask(0).shape()(0)==nx, AipsError);
152 : // AlwaysAssert(mask(0).shape()(1)==ny, AipsError);
153 0 : if((mask.shape()(0)!=nx) || (mask.shape()(1)!=ny)){
154 0 : throw(AipsError("Mask image shape is different from dirty image"));
155 : }
156 0 : if(nchan >1){
157 0 : if(mask.shape()(3)==nchan){
158 0 : isCubeMask=true;
159 0 : os << "Using multichannel mask" << LogIO::POST;
160 : }
161 : else{
162 : os << "Image cube and mask donot match in number of channels"
163 0 : << LogIO::WARN;
164 : os << "Will use first plane of the mask for all channels"
165 0 : << LogIO::WARN;
166 : }
167 : }
168 0 : LatticeStepper mls(mask.shape(),
169 0 : IPosition(4, nx, ny, 1, 1),
170 0 : IPosition(4, 0, 1, 3, 2));
171 0 : maskli= new RO_LatticeIterator<Float>(mask, mls);
172 0 : maskli->reset();
173 0 : mask_sl=makeMaskSubLat(nx, ny, newNx, newNy, *maskli, xbeg, xend,
174 : ybeg, yend);
175 0 : }
176 :
177 :
178 0 : maxresidual=0.0;
179 :
180 0 : for (Int ichan=0; ichan < nchan; ichan++) {
181 0 : if(useMask && isCubeMask && ichan >0) {
182 0 : (*maskli)++;
183 0 : if(mask_sl) delete mask_sl;
184 0 : mask_sl=0;
185 0 : mask_sl=makeMaskSubLat(nx, ny, newNx, newNy, *maskli, xbeg,
186 : xend, ybeg, yend);
187 : }
188 :
189 0 : if(nchan>1) {
190 0 : os<<"Processing channel "<<ichan+1<<" of "<<nchan<<LogIO::POST;
191 : }
192 0 : for (Int ipol=0; ipol < polloop; ++ipol){
193 0 : Int polbeg=0;
194 0 : Int polend=npol-1;
195 0 : if(!doPolJoint){
196 0 : polbeg=ipol;
197 0 : polend=ipol;
198 0 : os << "Doing stokes "<< stokesID(ipol) << " image" <<LogIO::POST;
199 : }
200 0 : LCBox imagebox(IPosition(4, xbeg, ybeg, polbeg, ichan),
201 0 : IPosition(4, xend, yend, polend, ichan),
202 0 : image.shape());
203 0 : LCBox psfbox(IPosition(4, 0, 0, 0, ichan),
204 0 : IPosition(4, nx-1, ny-1, 0, ichan),
205 0 : psf.shape());
206 :
207 0 : SubLattice<Float> psf_sl(psf, psfbox, false);
208 0 : SubLattice<Float> residual_sl(residual, imagebox, true);
209 0 : SubLattice<Float> model_sl=SubLattice<Float> (image, imagebox, true);
210 :
211 :
212 0 : ArrayLattice<Float> psftmp;
213 :
214 0 : if(nx != newNx){
215 :
216 0 : psftmp=ArrayLattice<Float> (IPosition(4, newNx, newNy, 1, 1));
217 0 : psftmp.set(0.0);
218 0 : Array<Float> tmparr=psf_sl.get();
219 0 : psftmp.putSlice(tmparr, IPosition(4, (newNx-nx)/2, (newNy-ny)/2, 0,0));
220 0 : psf_sl=SubLattice<Float>(psftmp, true);
221 :
222 0 : }
223 :
224 0 : TempLattice<Float> dirty_sl( residual_sl.shape());
225 0 : dirty_sl.copyData(residual_sl);
226 0 : TempLattice<Float> localmodel(model_sl.shape());
227 0 : localmodel.set(0.0);
228 :
229 : Float psfmax;
230 : {
231 0 : LatticeExprNode node = max(psf_sl);
232 0 : psfmax = node.getFloat();
233 0 : }
234 :
235 0 : if((psfmax==0.0) ||(useMask && (mask_sl == 0)) ) {
236 0 : maxresidual=0.0;
237 0 : iterused=0;
238 : os << LogIO::NORMAL // Loglevel INFO
239 0 : << "No data or blank mask for this channel: skipping" << LogIO::POST;
240 : } else {
241 0 : LatConvEquation eqn(psf_sl, residual_sl);
242 0 : ClarkCleanLatModel cleaner( localmodel );
243 0 : cleaner.setResidual(dirty_sl);
244 0 : if (mask_sl != 0 ) cleaner.setMask( *mask_sl );
245 :
246 : /*
247 : ClarkCleanProgress *cpp = 0;
248 : if (displayProgress_p) {
249 : cpp = new ClarkCleanProgress( pgplotter_p );
250 : cleaner.setProgress(*cpp);
251 : }
252 : */
253 :
254 0 : cleaner.setGain(gain);
255 0 : cleaner.setNumberIterations(numIter);
256 0 : cleaner.setThreshold(thresh);
257 0 : cleaner.setPsfPatchSize(IPosition(2,51,51));
258 0 : cleaner.setHistLength(1024);
259 0 : cleaner.setMaxNumPix(32*1024);
260 0 : cleaner.setCycleFactor(cycleFactor);
261 : // clean if there is no mask or if it has mask AND mask is not empty
262 0 : cleaner.solve(eqn);
263 0 : cleaner.setChoose(false);
264 0 : if (cleaner.numberIterations()>0) {
265 : os << LogIO::NORMAL // Loglevel INFO
266 : << "Clean used " << cleaner.numberIterations() << " iterations"
267 0 : << " in this round to get to a max residual of " << cleaner.threshold()
268 0 : << LogIO::POST;
269 : }
270 0 : LatticeExpr<Float> expr= model_sl + localmodel;
271 0 : model_sl.copyData(expr);
272 0 : maxresidual=max(maxresidual, cleaner.getMaxResidual());
273 0 : iterused=cleaner.numberIterations();
274 0 : converged = (cleaner.getMaxResidual() < thresh)
275 0 : || (cleaner.numberIterations()==0);
276 : // if (cpp != 0 ) delete cpp; cpp=0;
277 : // if (pgp != 0 ) delete pgp; pgp=0;
278 0 : }
279 0 : }
280 : }
281 0 : if (mask_sl != 0) {
282 0 : delete mask_sl;
283 0 : mask_sl=0;
284 : }
285 0 : if (maskli !=0) {
286 0 : delete maskli;
287 0 : maskli=0;
288 : }
289 : os << LogIO::NORMAL // Loglevel INFO
290 0 : << LatticeExprNode(sum(image)).getFloat()
291 : << " Jy <- The sum of the clean components"
292 0 : << LogIO::POST;
293 0 : return converged;
294 :
295 :
296 0 : }
297 :
298 0 : Lattice<Float>* ClarkCleanImageSkyModel::makeMaskSubLat(const Int& nx,
299 : const Int& ny,
300 : Int& newNx, Int& newNy,
301 : RO_LatticeIterator<Float>& maskIter,
302 : Int& xbeg,
303 : Int& xend,
304 : Int& ybeg,
305 : Int& yend) {
306 :
307 0 : LogIO os(LogOrigin("ClarkCleanImageSkyModel","makeMaskSubLat",WHERE));
308 :
309 :
310 0 : newNx=nx;
311 0 : newNy=ny;
312 0 : SubLattice<Float>* mask_sl = 0;
313 0 : xbeg=nx/4;
314 0 : ybeg=ny/4;
315 :
316 0 : xend=xbeg+nx/2-1;
317 0 : yend=ybeg+ny/2-1;
318 0 : Matrix<Float> mask= maskIter.matrixCursor();
319 0 : if(max(mask) < 0.000001) {
320 : //zero mask ....
321 0 : return mask_sl;
322 : }
323 : // Now read the mask and determine the bounding box
324 :
325 0 : xbeg=nx-1;
326 0 : ybeg=ny-1;
327 0 : xend=0;
328 0 : yend=0;
329 :
330 :
331 0 : for (Int iy=0;iy<ny;iy++) {
332 0 : for (Int ix=0;ix<nx;ix++) {
333 0 : if(mask(ix,iy)>0.000001) {
334 0 : xbeg=min(xbeg,ix);
335 0 : ybeg=min(ybeg,iy);
336 0 : xend=max(xend,ix);
337 0 : yend=max(yend,iy);
338 :
339 : }
340 : }
341 : }
342 : // Now have possible BLC. Make sure that we don't go over the
343 : // edge later
344 0 : Bool larger_quarter=false;
345 0 : if((xend - xbeg)>nx/2) {
346 : // xbeg=nx/4-1; //if larger than quarter take inner of mask
347 : //os << LogIO::WARN << "Mask span over more than half the x-axis: Considering inner half of the x-axis" << LogIO::POST;
348 0 : larger_quarter=true;
349 : }
350 0 : if((yend - ybeg)>ny/2) {
351 : //ybeg=ny/4-1;
352 : //os << LogIO::WARN << "Mask span over more than half the y-axis: Considering inner half of the y-axis" << LogIO::POST;
353 0 : larger_quarter=true;
354 : }
355 :
356 : // make sure xend, yend does not go out of bounds
357 0 : if(larger_quarter){
358 0 : xend=min(xend, nx-1);
359 0 : yend=min(yend, ny-1);
360 :
361 : }
362 : else{
363 0 : xend=min(xend,xbeg+nx/2-1);
364 0 : yend=min(yend,ybeg+ny/2-1);
365 : }
366 :
367 :
368 0 : if ((xend > xbeg) && (yend > ybeg) ) {
369 :
370 0 : IPosition latshape(4, mask.shape()(0), mask.shape()(1), 1,1);
371 0 : ArrayLattice<Float> arrayLat(latshape);
372 0 : LCBox maskbox (IPosition(4, xbeg, ybeg, 0, 0),
373 0 : IPosition(4, xend, yend, 0, 0),
374 0 : latshape);
375 :
376 0 : arrayLat.putSlice(mask, IPosition(4,0,0,0,0));
377 0 : mask_sl=new SubLattice<Float> (arrayLat, maskbox, false);
378 :
379 :
380 :
381 0 : }
382 :
383 0 : if(larger_quarter){
384 0 : newNx=2*nx;
385 0 : newNy=2*ny;
386 : /*xbeg=xbeg+(newNx-nx)/2;
387 : ybeg=ybeg+(newNy-ny)/2;
388 : xend=xend+(newNx-nx)/2;
389 : yend=yend+(newNy-ny)/2;
390 : */
391 : }
392 0 : return mask_sl;
393 0 : }
394 :
395 :
396 :
397 :
398 :
399 : } //# NAMESPACE CASA - END
400 :
|
