Line data Source code
1 : //# HogbomCleanImageSkyModel.cc: Implementation of HogbomCleanImageSkyModel 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/Cube.h>
29 : #include <casacore/casa/Arrays/Matrix.h>
30 : #include <casacore/casa/Arrays/ArrayMath.h>
31 : #include <casacore/casa/Arrays/ArrayMath.h>
32 : #include <casacore/casa/Arrays/ArrayPosIter.h>
33 : #include <msvis/MSVis/StokesVector.h>
34 : #include <synthesis/MeasurementComponents/HogbomCleanImageSkyModel.h>
35 : #include <casacore/images/Images/PagedImage.h>
36 : #include <casacore/casa/OS/File.h>
37 : #include <casacore/lattices/Lattices/LatticeStepper.h>
38 : #include <casacore/lattices/Lattices/LatticeIterator.h>
39 : #include <casacore/lattices/LEL/LatticeExpr.h>
40 : #include <casacore/lattices/LEL/LatticeExprNode.h>
41 : #include <synthesis/MeasurementEquations/SkyEquation.h>
42 : #include <casacore/casa/Exceptions/Error.h>
43 : #include <casacore/casa/BasicSL/String.h>
44 : #include <casacore/casa/Utilities/Assert.h>
45 :
46 : #include <sstream>
47 :
48 : #include <casacore/casa/Logging/LogMessage.h>
49 : #include <casacore/casa/Logging/LogSink.h>
50 :
51 : #include <casacore/casa/System/Choice.h>
52 :
53 :
54 : using namespace casacore;
55 : namespace casa { //# NAMESPACE CASA - BEGIN
56 :
57 : #define NEED_UNDERSCORES
58 : #if defined(NEED_UNDERSCORES)
59 : #define hclean hclean_
60 : #endif
61 : extern "C" {
62 : void hclean(Float*, Float*, Float*, int*, Float*, int*, int*, int*,
63 : int*, int*, int*, int*, int*, int*, int*, Float*, Float*,
64 : Float*, void *, void *);
65 : };
66 :
67 0 : void HogbomCleanImageSkyModelstopnow (Int *yes) {
68 0 : Vector<String> choices(2);
69 0 : choices(0)="Continue";
70 0 : choices(1)="Stop Now";
71 0 : LogMessage message(LogOrigin("HogbomCleanImageSkyModel","solve"));
72 0 : LogSink logSink;
73 0 : *yes=0;
74 0 : return;
75 : String choice=Choice::choice("Clean iteration: do you want to continue or stop?", choices);
76 : if (choice==choices(0)) {
77 : *yes=0;
78 : }
79 : else {
80 : message.message("Stopping");
81 : logSink.post(message);
82 : *yes=1;
83 : }
84 0 : }
85 :
86 0 : void HogbomCleanImageSkyModelmsgput(Int *npol, Int* /*pol*/, Int* iter, Int* px, Int* py,
87 : Float* fMaxVal) {
88 0 : LogMessage message(LogOrigin("HogbomCleanImageSkyModel","solve"));
89 0 : ostringstream o;
90 0 : LogSink logSink;
91 : /*
92 : String stokes("Unknown");
93 : ///UUU
94 : // This code allows only I, IV, IQU, IQUV !!!!!
95 : if(*npol==1) {
96 : stokes="I";
97 : }
98 : else if(*npol==2) {
99 : if(*pol==1) {
100 : stokes="I";
101 : }
102 : else if(*pol==2) {
103 : stokes="V";
104 : }
105 : }
106 : else if(*npol==3) {
107 : if(*pol==1) {
108 : stokes="I";
109 : }
110 : else if(*pol==2) {
111 : stokes="Q";
112 : }
113 : else if(*pol==3) {
114 : stokes="U";
115 : }
116 : }
117 : else if(*npol==4) {
118 : if(*pol==1) {
119 : stokes="I";
120 : }
121 : else if(*pol==2) {
122 : stokes="Q";
123 : }
124 : else if(*pol==3) {
125 : stokes="U";
126 : }
127 : else if(*pol==4) {
128 : stokes="V";
129 : }
130 : }
131 : if(*pol==-1) {
132 : stokes="I";
133 : }
134 : else if(*pol==-2) {
135 : stokes="I,V";
136 : }
137 : else if(*pol==-3) {
138 : stokes="I,Q,U";
139 : }
140 : else if(*pol==-4) {
141 : stokes="I,Q,U,V";
142 : }
143 : */
144 0 : if(*npol<0) {
145 0 : StokesVector maxVal(fMaxVal[0], fMaxVal[1], fMaxVal[2], fMaxVal[3]);
146 0 : if(*iter==0) {
147 : // o<<stokes<<": Before iteration, peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
148 0 : o<<"Before iteration, peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
149 0 : message.message(o);
150 0 : logSink.post(message);
151 : }
152 0 : else if(*iter>-1) {
153 : // o<<stokes<<": Iteration "<<*iter<<" peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
154 0 : o<<"Iteration "<<*iter<<" peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
155 0 : message.message(o);
156 0 : logSink.post(message);
157 : }
158 : else {
159 : // o<<stokes<<": Final iteration "<<abs(*iter)<<" peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
160 0 : o<<"Final iteration "<<abs(*iter)<<" peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
161 0 : message.message(o);
162 0 : logSink.post(message);
163 : }
164 : }
165 : else {
166 0 : Float maxVal(fMaxVal[0]);
167 0 : if(*iter==0) {
168 : // o<<stokes<<": Before iteration, peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
169 0 : o<<"Before iteration, peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
170 0 : message.message(o);
171 0 : logSink.post(message);
172 : }
173 0 : else if(*iter>-1) {
174 : // o<<stokes<<": Iteration "<<*iter<<" peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
175 0 : o<<"Iteration "<<*iter<<" peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
176 0 : message.message(o);
177 0 : logSink.post(message);
178 : }
179 : else {
180 : // o<<stokes<<": Final iteration "<<abs(*iter)<<" peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
181 0 : o<<"Final iteration "<<abs(*iter)<<" peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
182 0 : message.message(o);
183 0 : logSink.post(message);
184 : }
185 : }
186 :
187 0 : }
188 :
189 : // Clean solver
190 0 : Bool HogbomCleanImageSkyModel::solve(SkyEquation& se) {
191 :
192 0 : LogIO os(LogOrigin("HogbomCleanImageSkyModel","solve"));
193 :
194 0 : Bool converged=false;
195 0 : if(numberOfModels()>1) {
196 0 : os << "Cannot process more than one field" << LogIO::EXCEPTION;
197 : }
198 :
199 : // Make the residual image
200 0 : if(modified_p)
201 0 : makeNewtonRaphsonStep(se, false, (numberIterations()<1)?true:False);
202 :
203 : //Make the PSF
204 0 : if(!donePSF_p)
205 0 : makeApproxPSFs(se);
206 :
207 0 : if(numberIterations() < 1){
208 0 : return true;
209 : }
210 :
211 0 : if(!isSolveable(0)) {
212 0 : os << "Model 1 is not solveable!" << LogIO::EXCEPTION;
213 : }
214 :
215 :
216 :
217 0 : Int nx=image(0).shape()(0);
218 0 : Int ny=image(0).shape()(1);
219 0 : Int npol=image(0).shape()(2);
220 0 : Int nchan=image(0).shape()(3);
221 0 : Bool isCubeMask=false;
222 0 : AlwaysAssert((npol==1)||(npol==2)||(npol==3)||(npol==4), AipsError);
223 :
224 : // Loop over all channels
225 0 : LatticeStepper psfls(PSF(0).shape(), IPosition(4,nx,ny,1,1),
226 0 : IPosition(4,0,1,2,3));
227 0 : RO_LatticeIterator<Float> psfli(PSF(0),psfls);
228 :
229 : // Read the entire image for each spectral channel
230 0 : LatticeStepper ls(image(0).shape(),
231 0 : IPosition(4, nx, ny, npol, 1),
232 0 : IPosition(4, 0, 1, 2, 3));
233 0 : LatticeIterator<Float> imageStepli(residual(0), ls);
234 0 : LatticeIterator<Float> imageli(image(0), ls);
235 :
236 : // Now set up the mask.
237 0 : RO_LatticeIterator<Float>* maskli = 0;
238 0 : Matrix<Float>* lmask= new Matrix<Float>(1,1);
239 0 : lmask->set(1.0);
240 : Int xbeg, xend;
241 : Int ybeg, yend;
242 0 : Int newNx=nx;
243 0 : Int newNy=ny;
244 : //default clean support
245 0 : xbeg=nx/4;
246 0 : xend=3*nx/4-1;
247 0 : ybeg=ny/4;
248 0 : yend=3*ny/4-1;
249 0 : Bool domask=true;
250 0 : if(hasMask(0)) {
251 0 : domask=true;
252 0 : AlwaysAssert(mask(0).shape()(0)==nx, AipsError);
253 0 : AlwaysAssert(mask(0).shape()(1)==ny, AipsError);
254 0 : if(nchan >1){
255 0 : if(mask(0).shape()(3)==nchan){
256 0 : isCubeMask=true;
257 0 : os << "Using multichannel mask" << LogIO::POST;
258 : }
259 : else{
260 : os << "Image cube and mask donot match in number of channels"
261 0 : << LogIO::WARN;
262 : os << "Will use first plane of the mask for all channels"
263 0 : << LogIO::WARN;
264 : }
265 : }
266 0 : LatticeStepper mls(mask(0).shape(),
267 0 : IPosition(4, nx, ny, 1, 1),
268 0 : IPosition(4, 0, 1, 2, 3));
269 0 : maskli= new RO_LatticeIterator<Float>(mask(0), mls);
270 0 : maskli->reset();
271 0 : lmask=makeMaskMatrix(nx, ny, newNx, newNy, *maskli, xbeg,
272 : xend, ybeg, yend);
273 0 : }
274 : else {
275 0 : domask=false;
276 : }
277 :
278 : //Some variables if we change the image size
279 0 : IPosition newshp(4, newNx, newNy, npol,1);
280 0 : Array<Float> newData;
281 0 : Array<Float> newStep;
282 0 : Array<Float> newPsf;
283 : //------------------
284 :
285 0 : Float maxRes=0.0;
286 0 : Int iterUsed=0;
287 0 : Int chan=0;
288 0 : for (imageStepli.reset(),imageli.reset(),psfli.reset();
289 0 : !imageStepli.atEnd();
290 0 : imageStepli++,imageli++,psfli++,chan++) {
291 : //Deal with cube mask
292 0 : if(hasMask(0) && isCubeMask && chan >0) {
293 0 : (*maskli)++;
294 0 : if(lmask) delete lmask;
295 0 : lmask=0;
296 0 : lmask=makeMaskMatrix(nx, ny, newNx, newNy, *maskli, xbeg,
297 : xend, ybeg, yend);
298 : }
299 :
300 : // Make IPositions and find position of peak of PSF
301 0 : IPosition psfposmax(psfli.cursor().ndim());
302 0 : IPosition psfposmin(psfli.cursor().ndim());
303 : Float psfmax;
304 : Float psfmin;
305 0 : minMax(psfmin, psfmax, psfposmin, psfposmax, psfli.cursor());
306 :
307 0 : if(nchan>1) {
308 0 : os<<"Processing channel "<<chan+1<<" of "<<nchan<<LogIO::POST;
309 : }
310 0 : if((psfmax==0.0) || (hasMask(0) && (lmask==0))) {
311 0 : os<<"No data or blank mask for this channel: skipping"<<LogIO::POST;
312 : }
313 : else {
314 :
315 : Bool delete_iti, delete_its, delete_itp, delete_itm;
316 : const Float *lpsf_data, *lmask_data;
317 : Float *limage_data, *limageStep_data;
318 0 : lmask_data=lmask->getStorage(delete_itm);
319 0 : if(newNx==nx){
320 0 : limage_data=imageli.rwCursor().getStorage(delete_iti);
321 0 : limageStep_data=imageStepli.rwCursor().getStorage(delete_its);
322 0 : lpsf_data=psfli.cursor().getStorage(delete_itp);
323 : }
324 : else{
325 0 : IPosition blc(4, (newNx-nx)/2, (newNy-ny)/2, 0,0);
326 0 : IPosition trc(4, (newNx+nx)/2-1, (newNy+ny)/2-1, npol-1,0);
327 0 : newshp(0)=newNx;
328 0 : newshp(1)=newNy;
329 0 : newData.resize(newshp);
330 0 : newStep.resize(newshp);
331 :
332 0 : newData.set(0.0); newStep.set(0.0);
333 0 : newData(blc,trc).assign(imageli.rwCursor());
334 0 : newStep(blc,trc).assign(imageStepli.rwCursor());
335 : //psf we use first plane only of pol image
336 0 : newshp(2)=1;
337 0 : newPsf.resize(newshp);
338 0 : newPsf.set(0.0);
339 0 : trc(2)=0;
340 0 : newPsf(blc,trc).assign(psfli.cursor());
341 0 : limage_data=newData.getStorage(delete_iti);
342 0 : limageStep_data=newStep.getStorage(delete_its);
343 0 : lpsf_data=newPsf.getStorage(delete_itp);
344 0 : }
345 0 : Int niter=numberIterations();
346 0 : Float g=gain();
347 0 : Float thres=threshold();
348 : // Fortran indexes at 1
349 0 : Int fxbeg=xbeg+1;
350 0 : Int fxend=xend+1;
351 0 : Int fybeg=ybeg+1;
352 0 : Int fyend=yend+1;
353 : Int domaskI;
354 0 : if (domask) {
355 0 : domaskI = 1;
356 : } else {
357 0 : domaskI = 0;
358 : }
359 :
360 0 : Int starting_iteration = 0;
361 0 : Int ending_iteration=0;
362 : //# The const of lpsf and mask has to be cast away for the
363 : //# fortran function hclean.
364 0 : Float cycleSpeedup = -1; // ie, ignore it
365 0 : hclean(limage_data, limageStep_data,
366 : (Float*)lpsf_data, &domaskI, (Float*)lmask_data,
367 : &newNx, &newNy, &npol,
368 : &fxbeg, &fxend, &fybeg, &fyend, &niter,
369 : &starting_iteration, &ending_iteration,
370 : &g, &thres, &cycleSpeedup,
371 : (void*) &HogbomCleanImageSkyModelmsgput,
372 : (void*) &HogbomCleanImageSkyModelstopnow);
373 :
374 0 : iterUsed=max(iterUsed, ending_iteration);
375 0 : if(nx==newNx){
376 0 : imageli.rwCursor().putStorage (limage_data, delete_iti);
377 0 : imageStepli.rwCursor().putStorage (limageStep_data, delete_its);
378 0 : psfli.cursor().freeStorage (lpsf_data, delete_itp);
379 : }
380 : else{
381 0 : IPosition blc(4, (newNx-nx)/2, (newNy-ny)/2, 0,0);
382 0 : IPosition trc(4, (newNx+nx)/2-1, (newNy+ny)/2-1, npol-1,0);
383 0 : newData.putStorage(limage_data, delete_iti);
384 0 : newStep.putStorage(limageStep_data, delete_its);
385 0 : newPsf.freeStorage(lpsf_data, delete_itp);
386 0 : imageli.rwCursor().assign(newData(blc,trc));
387 0 : imageStepli.rwCursor().assign(newStep(blc,trc));
388 0 : }
389 0 : lmask->freeStorage (lmask_data, delete_itm);
390 : Float residualmax, residualmin;
391 0 : minMax(residualmin, residualmax, imageStepli.cursor());
392 0 : residualmax=max(residualmax, abs(residualmin));
393 0 : maxRes=max(maxRes, residualmax);
394 :
395 0 : converged = (residualmax < threshold());
396 : }
397 0 : if (lmask != 0 && isCubeMask) {
398 0 : lmask->resize(1,1);
399 : }
400 :
401 0 : }
402 0 : os << LatticeExprNode(sum(image(0))).getFloat()
403 0 : << " Jy is the sum of clean components " << LogIO::POST;
404 0 : modified_p=true;
405 0 : setThreshold(maxRes);
406 0 : setNumberIterations(iterUsed);
407 0 : if (maskli != 0) {
408 0 : delete maskli;
409 0 : maskli = 0;
410 : }
411 0 : return(converged);
412 0 : };
413 :
414 0 : Matrix<Float>* HogbomCleanImageSkyModel::makeMaskMatrix(const Int& nx,
415 : const Int& ny,
416 : Int& newNx, Int& newNy,
417 : RO_LatticeIterator<Float>& maskIter,
418 : Int& xbeg,
419 : Int& xend,
420 : Int& ybeg,
421 : Int& yend) {
422 :
423 0 : LogIO os(LogOrigin("HogbomCleanImageSkyModel","makeMaskMatrix",WHERE));
424 :
425 :
426 0 : newNx=nx;
427 0 : newNy=ny;
428 0 : xbeg=nx/4;
429 0 : ybeg=ny/4;
430 :
431 0 : xend=xbeg+nx/2-1;
432 0 : yend=ybeg+ny/2-1;
433 0 : Matrix<Float>* mask= new Matrix<Float>(maskIter.matrixCursor().shape());
434 0 : (*mask)=maskIter.matrixCursor();
435 : // ignore mask if none exists
436 0 : if(max(*mask) < 0.000001) {
437 : //os << "Mask seems to be empty; will CLEAN inner quarter"
438 : // << LogIO::WARN;
439 : //mask->resize(1,1);
440 : //mask->set(1.0);
441 0 : delete mask;
442 0 : mask=0;
443 0 : return mask;
444 : }
445 : // Now read the mask and determine the bounding box
446 0 : xbeg=nx-1;
447 0 : ybeg=ny-1;
448 0 : xend=0;
449 0 : yend=0;
450 :
451 0 : for (Int iy=0;iy<ny;iy++) {
452 0 : for (Int ix=0;ix<nx;ix++) {
453 0 : if((*mask)(ix,iy)>0.000001) {
454 0 : xbeg=min(xbeg,ix);
455 0 : ybeg=min(ybeg,iy);
456 0 : xend=max(xend,ix);
457 0 : yend=max(yend,iy);
458 :
459 : }
460 : }
461 : }
462 : // Now have possible BLC. Make sure that we don't go over the
463 : // edge later
464 0 : Bool larger_quarter=false;
465 0 : if((xend - xbeg)>nx/2) {
466 : //xbeg=nx/4-1; //if larger than quarter take inner of mask
467 : //os << LogIO::WARN << "Mask span over more than half the x-axis: Considering inner half of the x-axis" << LogIO::POST;
468 0 : larger_quarter=true;
469 : }
470 0 : if((yend - ybeg)>ny/2) {
471 : //ybeg=ny/4-1;
472 : //os << LogIO::WARN << "Mask span over more than half the y-axis: Considering inner half of the y-axis" << LogIO::POST;
473 0 : larger_quarter=true;
474 : }
475 0 : if(larger_quarter){
476 0 : newNx=2*nx;
477 0 : newNy=2*ny;
478 0 : Matrix<Float> * newMask= new Matrix<Float>(newNx, newNy);
479 0 : newMask->set(0.0);
480 0 : IPosition blc(2, (newNx-nx)/2, (newNy-ny)/2);
481 0 : IPosition trc(2, (newNx+nx)/2-1, (newNy+ny)/2-1);
482 0 : ((*newMask)(blc,trc)).assign(*mask);
483 0 : delete mask;
484 0 : mask=newMask;
485 0 : xbeg=xbeg+(newNx-nx)/2;
486 0 : ybeg=ybeg+(newNy-ny)/2;
487 0 : xend=xend+(newNx-nx)/2;
488 0 : yend=yend+(newNy-ny)/2;
489 0 : }
490 :
491 0 : xend=min(xend,xbeg+newNx/2-1);
492 0 : yend=min(yend,ybeg+newNy/2-1);
493 :
494 :
495 0 : return mask;
496 0 : }
497 :
498 :
499 : } //# NAMESPACE CASA - END
500 :
|
