Line data Source code
1 : //# SDAlgorithmClarkClean.cc: Implementation of SDAlgorithmClarkClean 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 :
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 : #include <synthesis/ImagerObjects/SDAlgorithmClarkClean.h>
61 : //#include <synthesis/MeasurementEquations/ClarkCleanLatModel.h>
62 : //#include <synthesis/MeasurementEquations/LatConvEquation.h>
63 :
64 : #include <synthesis/MeasurementEquations/ClarkCleanModel.h>
65 : #include <synthesis/MeasurementEquations/ConvolutionEquation.h>
66 :
67 : using namespace casacore;
68 : namespace casa { //# NAMESPACE CASA - BEGIN
69 :
70 0 : SDAlgorithmClarkClean::SDAlgorithmClarkClean(String clarktype):
71 0 : SDAlgorithmBase()
72 : {
73 0 : itsAlgorithmName=String(clarktype);
74 0 : itsMatDeltaModel.resize();
75 0 : }
76 :
77 0 : SDAlgorithmClarkClean::~SDAlgorithmClarkClean()
78 : {
79 :
80 0 : }
81 :
82 : // void SDAlgorithmClarkClean::initializeDeconvolver( Float &peakresidual, Float &modelflux )
83 0 : Long SDAlgorithmClarkClean::estimateRAM(const vector<int>& imsize){
84 0 : Long mem=0;
85 0 : if(itsImages)
86 : //Clark deconvolvers will have psf + residual + model + mask (1 plane at a time)
87 0 : mem=sizeof(Float)*(itsImages->getShape()(0))*(itsImages->getShape()(1))*6/1024;
88 0 : else if(imsize.size() >1)
89 0 : mem=sizeof(Float)*(imsize[0])*(imsize[1])*6/1024;
90 : else
91 0 : return 0;
92 : //throw(AipsError("Deconvolver cannot estimate the memory usage at this point"));
93 0 : return mem;
94 : }
95 0 : void SDAlgorithmClarkClean::initializeDeconvolver()
96 : {
97 0 : LogIO os( LogOrigin("SDAlgorithmClarkClean","initializeDeconvolver",WHERE) );
98 :
99 0 : itsImages->residual()->get( itsMatResidual, true );
100 0 : itsImages->model()->get( itsMatModel, true );
101 0 : itsImages->psf()->get( itsMatPsf, true );
102 0 : itsImages->mask()->get( itsMatMask, true );
103 :
104 : /*
105 : cout << "Clark : initDecon : " << itsImages->residual()->shape() << " : " << itsMatResidual.shape()
106 : << itsImages->model()->shape() << " : " << itsMatModel.shape()
107 : << itsImages->psf()->shape() << " : " << itsMatPsf.shape()
108 : << endl;
109 : */
110 : /*
111 : IPosition shp = itsMatResidual.shape();
112 : IPosition startp( shp.nelements(),0);
113 : IPosition stopp( shp.nelements(),0);
114 : stopp[0]=shp[0]-1; stopp[1]=shp[1]-1;
115 : Matrix<Float> oneplane;
116 : oneplane.doNotDegenerate( itsMatResidual(startp,stopp), IPosition( ) );
117 : findMaxAbs( oneplane, itsPeakResidual, itsMaxPos );
118 : */
119 :
120 : /*
121 : findMaxAbs( itsMatResidual, itsPeakResidual, itsMaxPos );
122 : itsModelFlux = sum( itsMatModel );
123 :
124 : peakresidual = itsPeakResidual;
125 : modelflux = itsModelFlux;
126 : */
127 :
128 : // Initialize the Delta Image model. Resize if needed.
129 0 : if ( itsMatDeltaModel.shape().nelements() != itsMatModel.shape().nelements() )
130 0 : { itsMatDeltaModel.resize ( itsMatModel.shape() ); }
131 :
132 :
133 : //cout << "Image Shapes : " << itsMatResidual.shape() << endl;
134 :
135 : ////////////////////////////////////////////////////////////////////////////////////
136 : // Get psf patch size
137 : // ...........fill this in....... after fitted beams are precomputed earlier..... and stored in ImageStore.
138 :
139 : // 4 pixels: pretty arbitrary, but only look for sidelobes
140 : // outside the inner (2n+1) * (2n+1) square
141 0 : Int ncent=4;
142 :
143 : {//locate peak size
144 0 : CoordinateSystem cs= itsImages->psf()->coordinates();
145 0 : Vector<String> unitas=cs.worldAxisUnits();
146 0 : unitas(0)="arcsec";
147 0 : unitas(1)="arcsec";
148 0 : cs.setWorldAxisUnits(unitas);
149 : //Get the minimum increment in arcsec
150 0 : Double incr=abs(min(cs.increment()(0), cs.increment()(1)));
151 0 : if(incr > 0.0){
152 0 : GaussianBeam beamModel=itsImages->getBeamSet()(0,0);
153 : // GaussianBeam beamModel=beam(model)(0,0);
154 0 : ncent=max(ncent, Int(ceil(beamModel.getMajor("arcsec")/incr)));
155 0 : ncent=max(ncent, Int(ceil(beamModel.getMinor("arcsec")/incr)));
156 : }
157 : }
158 :
159 0 : psfpatch_p=3*ncent+1;
160 :
161 0 : os << "Choosing a PSF patch size of " << psfpatch_p << " pixels."<< LogIO::POST;
162 :
163 0 : }
164 :
165 0 : void SDAlgorithmClarkClean::takeOneStep( Float loopgain,
166 : Int cycleNiter,
167 : Float cycleThreshold,
168 : Float &peakresidual,
169 : Float &modelflux,
170 : Int &iterdone)
171 : {
172 :
173 : //// Store current model in this matrix.
174 0 : itsImages->model()->get( itsMatDeltaModel, true );
175 0 : itsMatModel.assign( itsMatDeltaModel ); // This should make an explicit copy
176 :
177 : //// Set model to zero
178 0 : itsImages->model()->set( 0.0 );
179 :
180 : //// Prepare a single plane mask. For multi-stokes, it still needs a single mask.
181 0 : IPosition shp = itsImages->mask()->shape();
182 0 : IPosition startp( shp.nelements(),0);
183 0 : IPosition stopp( shp.nelements(),1);
184 0 : stopp[0]=shp[0]; stopp[1]=shp[1];
185 0 : Slicer oneslice(startp, stopp);
186 : //cout << "Making mask slice of : " << oneslice.start() << " : " << oneslice.end() << endl;
187 0 : SubImage<Float> oneplane( *(itsImages->mask()), oneslice ) ;
188 :
189 : //// Set up the cleaner
190 0 : ClarkCleanModel cleaner(itsMatDeltaModel);
191 0 : cleaner.setMask(oneplane.get());
192 0 : cleaner.setGain(loopgain);
193 0 : cleaner.setNumberIterations(cycleNiter);
194 0 : cleaner.setInitialNumberIterations(0);
195 0 : cleaner.setThreshold(cycleThreshold);
196 0 : cleaner.setPsfPatchSize(IPosition(2,psfpatch_p));
197 0 : cleaner.setMaxNumberMajorCycles(10);
198 0 : cleaner.setHistLength(1024);
199 0 : cleaner.setMaxNumPix(32*1024);
200 0 : cleaner.setChoose(false);
201 0 : cleaner.setCycleSpeedup(-1.0);
202 0 : cleaner.setSpeedup(0.0);
203 :
204 : /*
205 : cleaner.setMask( oneplane );
206 : cleaner.setNumberIterations(cycleNiter);
207 : cleaner.setMaxNumberMajorCycles(10);
208 : cleaner.setMaxNumberMinorIterations(cycleNiter);
209 : cleaner.setGain(loopgain);
210 : cleaner.setThreshold(cycleThreshold);
211 : cleaner.setPsfPatchSize(IPosition(2,psfpatch_p));
212 : cleaner.setHistLength(1024);
213 : // cleaner.setMaxExtPsf(..)
214 : cleaner.setSpeedup(-1.0);
215 : cleaner.setMaxNumPix(32*1024);
216 : */
217 :
218 : //// Prepare a single plane PSF
219 0 : IPosition shp2 = itsImages->psf()->shape();
220 0 : IPosition startp2( shp2.nelements(),0);
221 0 : IPosition stopp2( shp2.nelements(),1);
222 0 : stopp2[0]=shp2[0]; stopp2[1]=shp2[1];
223 0 : Slicer oneslice2(startp2, stopp2);
224 : //cout << "Making psf slice of : " << oneslice2.start() << " : " << oneslice2.end() << endl;
225 0 : SubImage<Float> oneplane2( *(itsImages->psf()), oneslice2 ) ;
226 :
227 : //// Make the convolution equation with a single plane PSF and possibly multiplane residual
228 0 : Array<Float> psfplane;
229 0 : oneplane2.get( psfplane, true );
230 0 : ConvolutionEquation eqn( psfplane , itsMatResidual );
231 :
232 : //Bool result =
233 0 : cleaner.singleSolve( eqn, itsMatResidual );
234 : // cout << result << endl;
235 :
236 0 : iterdone = cleaner.numberIterations();
237 :
238 : // Retrieve residual before major cycle
239 : // itsImages->residual()->copyData( cleaner.getResidual() );
240 : // (itsImages->residual())->put( itsMatResidual );
241 :
242 0 : cleaner.getModel( itsMatDeltaModel );
243 :
244 : // Add delta model to old model
245 : //Bool ret2 =
246 : // itsImages->model()->get( itsMatDeltaModel, true );
247 : // itsMatModel += itsMatDeltaModel;
248 :
249 0 : itsMatModel = itsMatDeltaModel;
250 :
251 : //--------------------------------- DECIDE WHICH PEAK RESIDUAL TO USE HERE.....
252 :
253 : ////// Find Peak Residual across the whole image
254 :
255 : //itsImages->residual()->get( itsMatResidual, true );
256 : // itsImages->mask()->get( itsMatMask, true );
257 : // findMaxAbsMask( itsMatResidual, itsMatMask, itsPeakResidual, itsMaxPos );
258 :
259 :
260 : ////// Find Peak Residual from the Clark Clean method (within current active pixels only).
261 0 : itsPeakResidual = cleaner.getMaxResidual();
262 :
263 0 : peakresidual = itsPeakResidual;
264 :
265 : // Find Total Model flux
266 0 : modelflux = sum( itsMatModel ); // Performance hog ?
267 0 : (itsImages->model())->put( itsMatModel );
268 0 : }
269 :
270 0 : void SDAlgorithmClarkClean::finalizeDeconvolver()
271 : {
272 0 : (itsImages->residual())->put( itsMatResidual );
273 0 : (itsImages->model())->put( itsMatModel );
274 0 : }
275 :
276 :
277 0 : void SDAlgorithmClarkClean::queryDesiredShape(Int &nchanchunks,
278 : Int &npolchunks,
279 : IPosition imshape) // , nImageFacets.
280 : {
281 0 : AlwaysAssert( imshape.nelements()==4, AipsError );
282 0 : nchanchunks=imshape(3); // Each channel should appear separately.
283 :
284 0 : itsAlgorithmName.downcase();
285 :
286 0 : if( itsAlgorithmName.matches("clarkstokes") )
287 : {
288 0 : npolchunks=imshape(2); // Each pol should appear separately.
289 : }
290 : else // "Clark"
291 : {
292 0 : npolchunks=1; // Send in all pols together.
293 : }
294 :
295 0 : }
296 :
297 :
298 : } //# NAMESPACE CASA - END
299 :
|
