Line data Source code
1 : //# SDAlgorithmHogbomClean.cc: Implementation of SDAlgorithmHogbomClean classes
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 : #include <casacore/casa/Arrays/ArrayMath.h>
29 : #include <casacore/casa/OS/HostInfo.h>
30 : #include <synthesis/ImagerObjects/SDAlgorithmHogbomClean.h>
31 : #include <components/ComponentModels/SkyComponent.h>
32 : #include <components/ComponentModels/ComponentList.h>
33 : #include <casacore/images/Images/TempImage.h>
34 : #include <casacore/images/Images/SubImage.h>
35 : #include <casacore/images/Regions/ImageRegion.h>
36 : #include <casacore/casa/OS/File.h>
37 : #include <casacore/lattices/LEL/LatticeExpr.h>
38 : #include <casacore/lattices/Lattices/TiledLineStepper.h>
39 : #include <casacore/lattices/Lattices/LatticeStepper.h>
40 : #include <casacore/lattices/Lattices/LatticeIterator.h>
41 : #include <synthesis/TransformMachines/StokesImageUtil.h>
42 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
43 : #include <casacore/casa/Exceptions/Error.h>
44 : #include <casacore/casa/BasicSL/String.h>
45 : #include <casacore/casa/Utilities/Assert.h>
46 : #include <casacore/casa/OS/Directory.h>
47 : #include <casacore/tables/Tables/TableLock.h>
48 :
49 : #include<synthesis/ImagerObjects/SIMinorCycleController.h>
50 :
51 : #include <sstream>
52 :
53 : #include <casacore/casa/Logging/LogMessage.h>
54 : #include <casacore/casa/Logging/LogIO.h>
55 : #include <casacore/casa/Logging/LogSink.h>
56 :
57 : #include <casacore/casa/System/Choice.h>
58 : #include <msvis/MSVis/StokesVector.h>
59 :
60 :
61 : using namespace casacore;
62 : namespace casa { //# NAMESPACE CASA - BEGIN
63 :
64 : #define NEED_UNDERSCORES
65 : #if defined(NEED_UNDERSCORES)
66 : #define hclean hclean_
67 : #endif
68 : extern "C" {
69 : void hclean(Float*, Float*, Float*, int*, Float*, int*, int*, int*,
70 : int*, int*, int*, int*, int*, int*, int*, Float*, Float*,
71 : Float*, void *, void *);
72 : };
73 :
74 :
75 :
76 : ////////////////// Global functions ?
77 0 : void REFHogbomCleanImageSkyModelstopnow (Int *yes) {
78 :
79 0 : *yes=0;
80 0 : return;
81 :
82 : /*
83 : Vector<String> choices(2);
84 : choices(0)="Continue";
85 : choices(1)="Stop Now";
86 : LogMessage message(LogOrigin("REFHogbomCleanImageSkyModel","solve"));
87 : LogSink logSink;
88 : *yes=0;
89 : return;
90 : String choice=Choice::choice("Clean iteration: do you want to continue or stop?", choices);
91 : if (choice==choices(0)) {
92 : *yes=0;
93 : }
94 : else {
95 : message.message("Stopping");
96 : logSink.post(message);
97 : *yes=1;
98 : }
99 : */
100 :
101 : }
102 :
103 0 : void REFHogbomCleanImageSkyModelmsgput(Int *npol, Int* /*pol*/, Int* iter, Int* px, Int* py,
104 : Float* fMaxVal) {
105 0 : LogIO os(LogOrigin("REFHogbomCleanImageSkyModel","solve"));
106 0 : ostringstream o;
107 : // LogSink logSink(LogMessage::NORMAL2);
108 :
109 0 : if(*npol<0) {
110 0 : StokesVector maxVal(fMaxVal[0], fMaxVal[1], fMaxVal[2], fMaxVal[3]);
111 0 : if(*iter==0) {
112 : // o<<stokes<<": Before iteration, peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
113 0 : o<<"Before iteration, peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
114 0 : os << LogIO::NORMAL1 << o.str() << LogIO::POST;
115 : // message.message(o);
116 : // logSink.post(message);
117 : }
118 0 : else if(*iter>-1) {
119 : // o<<stokes<<": Iteration "<<*iter<<" peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
120 0 : o<<"Iteration "<<*iter<<" peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
121 0 : os << LogIO::NORMAL1 << o.str() << LogIO::POST;
122 : // message.message(o);
123 : // logSink.post(message);
124 : }
125 : else {
126 : // o<<stokes<<": Final iteration "<<abs(*iter)<<" peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
127 0 : o<<"Final iteration "<<abs(*iter)<<" peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
128 0 : os << LogIO::NORMAL1 << o.str() << LogIO::POST;
129 : // message.message(o);
130 : // logSink.post(message);
131 : }
132 : }
133 : else {
134 0 : Float maxVal(fMaxVal[0]);
135 0 : if(*iter==0) {
136 : // o<<stokes<<": Before iteration, peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
137 0 : o<<"Before iteration, peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
138 0 : os << LogIO::NORMAL1 << o.str() << LogIO::POST;
139 : // message.message(o);
140 : // logSink.post(message);
141 : }
142 0 : else if(*iter>-1) {
143 : // o<<stokes<<": Iteration "<<*iter<<" peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
144 0 : o<<"Iteration "<<*iter<<" peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
145 0 : os << LogIO::NORMAL1 << o.str() << LogIO::POST;
146 : // message.message(o);
147 : // logSink.post(message);
148 : }
149 : else {
150 : // o<<stokes<<": Final iteration "<<abs(*iter)<<" peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
151 0 : o<<"Final iteration "<<abs(*iter)<<" peak is "<<maxVal<<" at "<<*px-1<<","<<*py-1;
152 0 : os << LogIO::NORMAL1 << o.str() << LogIO::POST;
153 : // message.message(o);
154 : // logSink.post(message);
155 : }
156 : }
157 :
158 0 : }
159 :
160 :
161 :
162 0 : SDAlgorithmHogbomClean::SDAlgorithmHogbomClean():
163 0 : SDAlgorithmBase()
164 : // itsMatResidual(), itsMatModel(), itsMatPsf(),
165 : // itsMaxPos( IPosition() ),
166 : // itsPeakResidual(0.0),
167 : // itsModelFlux(0.0),
168 : // itsMatMask()
169 : {
170 0 : itsAlgorithmName=String("Hogbom");
171 0 : }
172 :
173 0 : SDAlgorithmHogbomClean::~SDAlgorithmHogbomClean()
174 : {
175 :
176 0 : }
177 :
178 : // void SDAlgorithmHogbomClean::initializeDeconvolver( Float &peakresidual, Float &modelflux )
179 0 : void SDAlgorithmHogbomClean::initializeDeconvolver()
180 : {
181 0 : LogIO os( LogOrigin("SDAlgorithmHogbomClean","initializeDeconvolver",WHERE) );
182 :
183 0 : itsImages->residual()->get( itsMatResidual, true );
184 0 : itsImages->model()->get( itsMatModel, true );
185 0 : itsImages->psf()->get( itsMatPsf, true );
186 0 : itsImages->mask()->get( itsMatMask, true );
187 :
188 : // cout << "initDecon : " << itsImages->residual()->shape() << " : " << itsMatResidual.shape()
189 : // << itsImages->model()->shape() << " : " << itsMatModel.shape()
190 : // << itsImages->psf()->shape() << " : " << itsMatPsf.shape()
191 : // << endl;
192 :
193 : /*
194 : findMaxAbs( itsMatResidual, itsPeakResidual, itsMaxPos );
195 : itsModelFlux = sum( itsMatModel );
196 :
197 : peakresidual = itsPeakResidual;
198 : modelflux = itsModelFlux;
199 : */
200 0 : }
201 :
202 :
203 0 : void SDAlgorithmHogbomClean::takeOneStep( Float loopgain,
204 : Int cycleNiter,
205 : Float cycleThreshold,
206 : Float &peakresidual,
207 : Float &modelflux,
208 : Int &iterdone)
209 : {
210 :
211 : Bool delete_iti, delete_its, delete_itp, delete_itm;
212 : const Float *lpsf_data, *lmask_data;
213 : Float *limage_data, *limageStep_data;
214 :
215 0 : limage_data = itsMatModel.getStorage( delete_iti );
216 0 : limageStep_data = itsMatResidual.getStorage( delete_its );
217 0 : lpsf_data = itsMatPsf.getStorage( delete_itp );
218 0 : lmask_data = itsMatMask.getStorage( delete_itm );
219 :
220 0 : Int niter= cycleNiter;
221 0 : Float g = loopgain;
222 0 : Float thres = cycleThreshold;
223 :
224 0 : IPosition shp = itsMatPsf.shape();
225 : /*
226 : Int xbeg = shp[0]/4;
227 : Int xend = 3*shp[0]/4;
228 : Int ybeg = shp[1]/4;
229 : Int yend = 3*shp[1]/4;
230 : */
231 :
232 0 : Int xbeg = 0;
233 0 : Int xend = shp[0]-1;
234 0 : Int ybeg = 0;
235 0 : Int yend = shp[1]-1;
236 :
237 :
238 0 : Int newNx = shp[0];
239 0 : Int newNy = shp[1];
240 0 : Int npol = 1;
241 :
242 : // Fortran indexes at 1
243 0 : Int fxbeg=xbeg+1;
244 0 : Int fxend=xend+1;
245 0 : Int fybeg=ybeg+1;
246 0 : Int fyend=yend+1;
247 :
248 0 : Int domaskI = 1;
249 :
250 0 : Int starting_iteration = 0;
251 0 : Int ending_iteration=0;
252 0 : Float cycleSpeedup = -1; // ie, ignore it
253 :
254 0 : hclean(limage_data, limageStep_data,
255 : (Float*)lpsf_data, &domaskI, (Float*)lmask_data,
256 : &newNx, &newNy, &npol,
257 : &fxbeg, &fxend, &fybeg, &fyend, &niter,
258 : &starting_iteration, &ending_iteration,
259 : &g, &thres, &cycleSpeedup,
260 : (void*) &REFHogbomCleanImageSkyModelmsgput,
261 : (void*) &REFHogbomCleanImageSkyModelstopnow);
262 :
263 0 : iterdone=ending_iteration;
264 :
265 0 : itsMatModel.putStorage( limage_data, delete_iti );
266 0 : itsMatResidual.putStorage( limageStep_data, delete_its );
267 0 : itsMatPsf.freeStorage( lpsf_data, delete_itp );
268 0 : itsMatMask.freeStorage( lmask_data, delete_itm );
269 :
270 :
271 : /////////////////
272 0 : findMaxAbsMask( itsMatResidual, itsMatMask, itsPeakResidual, itsMaxPos );
273 0 : peakresidual = itsPeakResidual;
274 :
275 0 : modelflux = sum( itsMatModel ); // Performance hog ?
276 0 : }
277 :
278 0 : void SDAlgorithmHogbomClean::finalizeDeconvolver()
279 : {
280 0 : (itsImages->residual())->put( itsMatResidual );
281 0 : (itsImages->model())->put( itsMatModel );
282 0 : }
283 :
284 :
285 : } //# NAMESPACE CASA - END
286 :
|
