Line data Source code
1 : //# SDAlgorithmClarkClean2.cc: Implementation of SDAlgorithmClarkClean2 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: 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 <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/SDAlgorithmClarkClean2.h>
61 : #include <synthesis/MeasurementEquations/ClarkCleanLatModel.h>
62 : #include <synthesis/MeasurementEquations/LatConvEquation.h>
63 :
64 :
65 : using namespace casacore;
66 : namespace casa { //# NAMESPACE CASA - BEGIN
67 :
68 0 : SDAlgorithmClarkClean2::SDAlgorithmClarkClean2(String clarktype):
69 0 : SDAlgorithmBase()
70 : {
71 0 : itsAlgorithmName=String(clarktype);
72 0 : itsMatDeltaModel.resize();
73 0 : itsMatModel.resize();
74 0 : }
75 :
76 0 : SDAlgorithmClarkClean2::~SDAlgorithmClarkClean2()
77 : {
78 :
79 0 : }
80 0 : Long SDAlgorithmClarkClean2::estimateRAM(const vector<int>& imsize){
81 0 : Long mem=0;
82 0 : if(itsImages)
83 : //Clark deconvolvers will have psf + residual + model + mask (1 plane at a time)
84 0 : mem=sizeof(Float)*(itsImages->getShape()(0))*(itsImages->getShape()(1))*6/1024;
85 0 : else if(imsize.size() >1)
86 0 : mem=sizeof(Float)*(imsize[0])*(imsize[1])*6/1024;
87 : else
88 0 : return 0;
89 : //throw(AipsError("Deconvolver cannot estimate the memory usage at this point"));
90 0 : return mem;
91 : }
92 :
93 : // void SDAlgorithmClarkClean2::initializeDeconvolver( Float &peakresidual, Float &modelflux )
94 0 : void SDAlgorithmClarkClean2::initializeDeconvolver()
95 : {
96 0 : LogIO os( LogOrigin("SDAlgorithmClarkClean2","initializeDeconvolver",WHERE) );
97 :
98 0 : itsImages->residual()->get( itsMatResidual, true );
99 :
100 : //itsImages->model()->get( itsMatModel, true ); not necessary here as it is done in Onestep all over again
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 : // if ( itsMatDeltaModel.shape().nelements() != itsMatModel.shape().nelements() )
130 : // { 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 0 : }
157 0 : }
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 SDAlgorithmClarkClean2::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 : LatticeLocker lockmod (*(itsImages->model()), FileLocker::Write);
175 0 : LatticeLocker lockresid (*(itsImages->residual()), FileLocker::Write);
176 0 : itsImages->model()->get( itsMatDeltaModel, true );
177 : // cerr << "Orig deltaMod " << max(itsMatDeltaModel) << endl;
178 0 : itsMatModel.resize();
179 0 : itsMatModel.assign( itsMatDeltaModel ); // This should make an explicit copy
180 : //cerr << "Orig Mod " << max(itsMatModel) << endl;
181 : //// Set model to zero
182 0 : itsImages->model()->set( 0.0 );
183 :
184 : // cerr << "max after set " << max(itsImages->model()->get(true)) << endl;
185 :
186 : //// Prepare a single plane mask. For multi-stokes, it still needs a single mask.
187 0 : IPosition shp = itsImages->mask()->shape();
188 0 : IPosition startp( shp.nelements(),0);
189 0 : IPosition stopp( shp.nelements(),1);
190 0 : stopp[0]=shp[0]; stopp[1]=shp[1];
191 0 : Slicer oneslice(startp, stopp);
192 : //cout << "Making mask slice of : " << oneslice.start() << " : " << oneslice.end() << endl;
193 0 : SubImage<Float> oneplane( *(itsImages->mask()), oneslice ) ;
194 :
195 : //// Set up the cleaner
196 0 : ClarkCleanLatModel cleaner(*(itsImages->model()));
197 0 : cleaner.setMask( oneplane );
198 0 : cleaner.setNumberIterations(cycleNiter);
199 0 : cleaner.setMaxNumberMajorCycles(10);
200 0 : cleaner.setMaxNumberMinorIterations(cycleNiter);
201 0 : cleaner.setGain(loopgain);
202 0 : cleaner.setThreshold(cycleThreshold);
203 0 : cleaner.setPsfPatchSize(IPosition(2,psfpatch_p));
204 0 : cleaner.setHistLength(1024);
205 : // cleaner.setMaxExtPsf(..)
206 0 : cleaner.setSpeedup(-1.0);
207 0 : cleaner.setMaxNumPix(32*1024);
208 :
209 : //// Prepare a single plane PSF
210 0 : IPosition shp2 = itsImages->psf()->shape();
211 0 : IPosition startp2( shp2.nelements(),0);
212 0 : IPosition stopp2( shp2.nelements(),1);
213 0 : stopp2[0]=shp2[0]; stopp2[1]=shp2[1];
214 0 : Slicer oneslice2(startp2, stopp2);
215 : // cout << "Making psf slice of : " << oneslice2.start() << " : " << oneslice2.end() << endl;
216 0 : SubImage<Float> oneplane2( *(itsImages->psf()), oneslice2 ) ;
217 :
218 : //// Make the convolution equation with a single plane PSF and possibly multiplane residual
219 0 : LatConvEquation eqn( oneplane2 , *(itsImages->residual()) );
220 :
221 : //Bool result =
222 0 : cleaner.solve( eqn );
223 : // cout << result << endl;
224 :
225 0 : iterdone = cleaner.getNumberIterations();
226 :
227 : // Retrieve residual before major cycle
228 0 : itsImages->residual()->copyData( cleaner.getResidual() );
229 0 : itsImages->residual()->get( itsMatResidual, true );
230 :
231 : // Add delta model to old model
232 : //Bool ret2 =
233 0 : itsImages->model()->get( itsMatDeltaModel, true );
234 : //cerr << "deltaModel " << max(itsMatDeltaModel) << endl;
235 0 : itsMatModel += itsMatDeltaModel;
236 : //cerr << "totModel " << max(itsMatModel) << endl;
237 : //--------------------------------- DECIDE WHICH PEAK RESIDUAL TO USE HERE.....
238 :
239 : ////// Find Peak Residual across the whole image
240 : /*
241 : itsImages->residual()->get( itsMatResidual, true );
242 : itsImages->mask()->get( itsMatMask, true );
243 : findMaxAbsMask( itsMatResidual, itsMatMask, itsPeakResidual, itsMaxPos );
244 : */
245 :
246 : ////// Find Peak Residual from the Clark Clean method (within current active pixels only).
247 0 : itsPeakResidual = cleaner.getMaxResidual();
248 :
249 0 : peakresidual = itsPeakResidual;
250 :
251 : // Find Total Model flux
252 0 : modelflux = sum( itsMatModel ); // Performance hog ?
253 : //cerr << "tot flux " << modelflux << endl;
254 0 : (itsImages->model())->put( itsMatModel );
255 0 : }
256 :
257 0 : void SDAlgorithmClarkClean2::finalizeDeconvolver()
258 : {
259 0 : (itsImages->residual())->put( itsMatResidual );
260 0 : (itsImages->model())->put( itsMatModel );
261 0 : }
262 :
263 :
264 0 : void SDAlgorithmClarkClean2::queryDesiredShape(Int &nchanchunks,
265 : Int &npolchunks,
266 : IPosition imshape) // , nImageFacets.
267 : {
268 0 : AlwaysAssert( imshape.nelements()==4, AipsError );
269 0 : nchanchunks=imshape(3); // Each channel should appear separately.
270 :
271 0 : itsAlgorithmName.downcase();
272 :
273 0 : if( itsAlgorithmName.matches("clarkstokes") )
274 : {
275 0 : npolchunks=imshape(2); // Each pol should appear separately.
276 : }
277 : else // "Clark"
278 : {
279 0 : npolchunks=1; // Send in all pols together.
280 : }
281 :
282 0 : }
283 :
284 :
285 : } //# NAMESPACE CASA - END
286 :
|
