Line data Source code
1 : //# SDAlgorithmBase.cc: Implementation of SDAlgorithmBase 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 : #include <synthesis/ImagerObjects/SDAlgorithmBase.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 : #include <imageanalysis/ImageAnalysis/CasaImageBeamSet.h>
51 :
52 : #include <sstream>
53 : #include <casacore/casa/OS/EnvVar.h>
54 :
55 : #include <casacore/casa/Logging/LogMessage.h>
56 : #include <casacore/casa/Logging/LogIO.h>
57 : #include <casacore/casa/Logging/LogSink.h>
58 :
59 : using namespace casacore;
60 : namespace casa { //# NAMESPACE CASA - BEGIN
61 :
62 :
63 1063 : SDAlgorithmBase::SDAlgorithmBase():
64 1063 : itsAlgorithmName("Test")
65 : // itsDecSlices (),
66 : // itsResidual(), itsPsf(), itsModel()
67 : {
68 1063 : }
69 :
70 1063 : SDAlgorithmBase::~SDAlgorithmBase()
71 : {
72 :
73 1063 : }
74 :
75 :
76 934 : void SDAlgorithmBase::deconvolve( SIMinorCycleController &loopcontrols,
77 : std::shared_ptr<SIImageStore> &imagestore,
78 : Int deconvolverid,
79 : Bool isautomasking, Bool fastnoise, Record robuststats, bool fullsummary)
80 : {
81 1868 : LogIO os( LogOrigin("SDAlgorithmBase","deconvolve",WHERE) );
82 :
83 : // Make a list of Slicers.
84 : //partitionImages( imagestore );
85 :
86 : Int nSubChans, nSubPols;
87 934 : queryDesiredShape(nSubChans, nSubPols, imagestore->getShape());
88 :
89 : // cout << "Check imstore from deconvolver : " << endl;
90 934 : imagestore->validate();
91 :
92 : //itsImages = imagestore;
93 :
94 : // loadMask();
95 :
96 : //os << "-------------------------------------------------------------------------------------------------------------" << LogIO::POST;
97 :
98 :
99 934 : os << "[" << imagestore->getName() << "]";
100 934 : os << " Run " << itsAlgorithmName << " minor-cycle " ;
101 934 : if( nSubChans>1 || nSubPols>1 ) os << "on ";
102 934 : if( nSubChans >1 ) os << nSubChans << " chans " ;
103 934 : if( nSubPols >1 ) os << nSubPols << " pols ";
104 1868 : os << "| CycleThreshold=" << loopcontrols.getCycleThreshold()
105 : << ", CycleNiter=" << loopcontrols.getCycleNiter()
106 1868 : << ", Gain=" << loopcontrols.getLoopGain()
107 934 : << LogIO::POST;
108 :
109 934 : Float itsPBMask = loopcontrols.getPBMask();
110 934 : Int cycleStartIteration = loopcontrols.getIterDone();
111 :
112 934 : Float maxResidualAcrossPlanes=0.0; Int maxResChan=0,maxResPol=0;
113 934 : Float totalFluxAcrossPlanes=0.0;
114 :
115 3870 : for( Int chanid=0; chanid<nSubChans;chanid++)
116 : {
117 6095 : for( Int polid=0; polid<nSubPols; polid++)
118 : {
119 : // itsImages = imagestore->getSubImageStoreOld( chanid, onechan, polid, onepol );
120 3159 : itsImages = imagestore->getSubImageStore( 0, 1, chanid, nSubChans, polid, nSubPols );
121 :
122 3159 : Int startiteration = loopcontrols.getIterDone(); // TODO : CAS-8767 key off subimage index
123 3159 : Float peakresidual=0.0, peakresidualnomask=0.0;
124 3159 : Float modelflux=0.0;
125 3159 : Int iterdone=0;
126 :
127 : ///itsMaskHandler.resetMask( itsImages ); //, (loopcontrols.getCycleThreshold()/peakresidual) );
128 3159 : Int stopCode=0;
129 :
130 3159 : Float startpeakresidual = 0.0, startpeakresidualnomask = 0.0;
131 3159 : Float startmodelflux = 0.0;
132 3159 : Array<Double> robustrms;
133 :
134 3159 : Float masksum = itsImages->getMaskSum();
135 3159 : Bool validMask = ( masksum > 0 );
136 :
137 3159 : peakresidualnomask = itsImages->getPeakResidual();
138 3159 : if( validMask ) peakresidual = itsImages->getPeakResidualWithinMask();
139 101 : else peakresidual = 0; // CAS-14201 : If mask is zero, peakresidual is zero
140 : //else peakresidual = peakresidualnomask;
141 3159 : modelflux = itsImages->getModelFlux();
142 :
143 3157 : startpeakresidual = peakresidual;
144 3157 : startpeakresidualnomask = peakresidualnomask;
145 3157 : startmodelflux = modelflux;
146 :
147 : //Float nsigma = 150.0; // will set by user, fixed for 3sigma for now.
148 3157 : Float nsigma = loopcontrols.getNsigma();
149 : //os<<"robustrms nelements="<<robustrms.nelements()<<LogIO::POST;
150 : Float nsigmathresh;
151 3157 : if (robustrms.nelements()==0) {
152 3157 : nsigmathresh = 0.0;
153 : } else{
154 0 : nsigmathresh = nsigma * (Float)robustrms(IPosition(1,0));
155 : }
156 :
157 : Float thresholdtouse;
158 3157 : if (nsigma>0.0) {
159 : // returns as an Array but itsImages is already single plane so
160 : // the return rms contains only a single element
161 385 : Array<Double> medians;
162 : //os<<"robuststats rec size="<<robuststats.nfields()<<LogIO::POST;
163 385 : Bool statsexists = false;
164 385 : IPosition statsindex;
165 385 : if (robuststats.nfields()) {
166 : // use existing stats
167 385 : if (robuststats.isDefined("robustrms")) {
168 385 : robuststats.get(RecordFieldId("robustrms"), robustrms);
169 385 : robuststats.get(RecordFieldId("median"), medians);
170 385 : statsexists=True;
171 :
172 : }
173 0 : else if(robuststats.isDefined("medabsdevmed")) {
174 0 : Array<Double> mads;
175 0 : robuststats.get(RecordFieldId("medabsdevmed"), mads);
176 0 : robuststats.get(RecordFieldId("median"), medians);
177 0 : robustrms = mads * 1.4826; // convert to rms
178 0 : statsexists=True;
179 :
180 0 : }
181 : //statsindex=IPosition(1,chanid); // this only support for npol =1, need to fix this
182 385 : if((robustrms.shape().nelements() ==2) && (robustrms.shape()[0] > polid) && (robustrms.shape()[1] > chanid) )
183 0 : statsindex=IPosition(2,polid,chanid);
184 385 : else if((robustrms.shape().nelements() ==1) && (robustrms.shape()[0] > chanid))
185 385 : statsindex=IPosition(1,chanid);
186 : else ///no idea what got passed in the record
187 0 : statsexists=False;
188 :
189 : }
190 385 : if (statsexists) {
191 385 : os<<LogIO::DEBUG1<<"Using the existing robust image statatistics!"<<LogIO::POST;
192 : }
193 : else {
194 0 : robustrms = itsImages->calcRobustRMS(medians, itsPBMask, fastnoise);
195 0 : statsindex=IPosition(1,0);
196 : }
197 385 : if (isautomasking) { // new threshold defination
198 : //nsigmathresh = (Float)medians(IPosition(1,0)) + nsigma * (Float)robustrms(IPosition(1,0));
199 4 : nsigmathresh = (Float)medians(statsindex) + nsigma * (Float)robustrms(statsindex);
200 : }
201 : else {
202 381 : nsigmathresh = nsigma * (Float)robustrms(statsindex);
203 : }
204 385 : thresholdtouse = max( nsigmathresh, loopcontrols.getCycleThreshold());
205 385 : }
206 : else {
207 2772 : thresholdtouse = loopcontrols.getCycleThreshold();
208 : }
209 : //os << LogIO::DEBUG1<<"loopcontrols.getCycleThreshold()="<<loopcontrols.getCycleThreshold()<<LogIO::POST;
210 3157 : os << LogIO::NORMAL3<<"current CycleThreshold="<<loopcontrols.getCycleThreshold()<<" nsigma threshold="<<nsigmathresh<<LogIO::POST;
211 3157 : String thresholddesc = (thresholdtouse == loopcontrols.getCycleThreshold() ? "cyclethreshold" : "n-sigma");
212 3157 : os << LogIO::NORMAL3<< "thresholdtouse="<< thresholdtouse << "("<<thresholddesc<<")"<< LogIO::POST;
213 :
214 3157 : if (thresholddesc=="n-sigma") {
215 : //os << LogIO::DEBUG1<< "Set nsigma thresh="<<nsigmathresh<<LogIO::POST;
216 31 : loopcontrols.setNsigmaThreshold(nsigmathresh);
217 : }
218 3157 : loopcontrols.setPeakResidual( peakresidual );
219 3157 : loopcontrols.resetMinResidual(); // Set it to current initial peakresidual.
220 :
221 3157 : stopCode = checkStop( loopcontrols, peakresidual );
222 :
223 3157 : if( validMask && stopCode==0 )
224 : {
225 : // Init the deconvolver
226 : //where we write in model and residual may be
227 : {
228 :
229 2837 : LatticeLocker lockresid (*(itsImages->residual()), FileLocker::Read);
230 2837 : LatticeLocker lockmodel (*(itsImages->model()), FileLocker::Read);
231 2837 : LatticeLocker lockmask (*(itsImages->mask()), FileLocker::Read);
232 2837 : LatticeLocker lockpsf (*(itsImages->psf()), FileLocker::Read);
233 :
234 2837 : initializeDeconvolver();
235 2855 : }
236 :
237 5662 : while ( stopCode==0 )
238 : {
239 :
240 2831 : if (nsigma>0.0) {
241 371 : os << "Using " << thresholddesc << " for threshold criterion: (cyclethreshold="<<loopcontrols.getCycleThreshold()<< ", nsigma threshold="<<nsigmathresh<<" )" << LogIO::POST;
242 371 : loopcontrols.setNsigmaThreshold(nsigmathresh);
243 : }
244 2831 : Int thisniter = loopcontrols.getCycleNiter() <5000 ? loopcontrols.getCycleNiter() : 2000;
245 :
246 2831 : loopcontrols.setPeakResidual( peakresidual );
247 :
248 : //where we write in model and residual may be
249 : {
250 : //no need to lock here as only arrays are used
251 :
252 : //LatticeLocker lock1 (*(itsImages->residual()), FileLocker::Write);
253 : //LatticeLocker lock2 (*(itsImages->model()), FileLocker::Write);
254 2831 : takeOneStep( loopcontrols.getLoopGain(),
255 : // loopcontrols.getCycleNiter(),
256 : thisniter,
257 : //loopcontrols.getCycleThreshold(),
258 : thresholdtouse,
259 : peakresidual,
260 : modelflux,
261 : iterdone);
262 : }
263 :
264 2831 : os << LogIO::NORMAL1 << "SDAlgoBase: After one step, dec : " << deconvolverid << " residual=" << peakresidual << " model=" << modelflux << " iters=" << iterdone << LogIO::POST;
265 :
266 2831 : SynthesisUtilMethods::getResource("In Deconvolver : one step" );
267 :
268 2831 : loopcontrols.incrementMinorCycleCount( iterdone ); // CAS-8767 : add subimageindex and merge with addSummaryMinor call later.
269 :
270 2831 : stopCode = checkStop( loopcontrols, peakresidual );
271 :
272 : /// Catch the situation where takeOneStep returns without satisfying any
273 : /// convergence criterion. For now, takeOneStep is the entire minor cycle.
274 : /// Later, when you can interrupt minor cycles, takeOneStep will become more
275 : /// fine grained, and then stopCode=0 will be valid. For now though, check on
276 : /// it and handle it (for CAS-7898).
277 2831 : if(stopCode==0 && iterdone != thisniter)
278 : {
279 7 : os << LogIO::NORMAL1 << "Exited " << itsAlgorithmName << " minor cycle without satisfying stopping criteria " << LogIO::POST;
280 7 : stopCode=5;
281 : }
282 :
283 : }// end of minor cycle iterations for this subimage.
284 :
285 : //where we write in model and residual may be
286 : {
287 :
288 2831 : LatticeLocker lock1 (*(itsImages->residual()), FileLocker::Write);
289 2831 : LatticeLocker lock2 (*(itsImages->model()), FileLocker::Write);
290 2831 : finalizeDeconvolver();
291 2831 : }
292 :
293 : // get the new peakres without a mask for the summary minor
294 2831 : peakresidualnomask = itsImages->getPeakResidual();
295 :
296 : }// if validmask
297 :
298 : // same as checking on itscycleniter.....
299 3151 : loopcontrols.setUpdatedModelFlag( loopcontrols.getIterDone()-startiteration );
300 :
301 3151 : os << "[" << imagestore->getName();
302 3151 : if(nSubChans>1) os << ":C" << chanid ;
303 3151 : if(nSubPols>1) os << ":P" << polid ;
304 3151 : Int iterend = loopcontrols.getIterDone();
305 : os << "]"
306 : // <<" iters=" << ( (iterend>startiteration) ? startiteration+1 : startiteration )<< "->" << iterend
307 : <<" iters=" << startiteration << "->" << iterend << " [" << iterend-startiteration << "]"
308 : << ", model=" << startmodelflux << "->" << modelflux
309 3151 : << ", peakres=" << startpeakresidual << "->" << peakresidual ;
310 :
311 3151 : switch (stopCode)
312 : {
313 0 : case 0:
314 0 : os << ", Skipped this plane. Zero mask.";
315 0 : break;
316 1911 : case 1:
317 1911 : os << ", Reached cycleniter.";
318 1911 : break;
319 1199 : case 2:
320 1199 : os << ", Reached cyclethreshold.";
321 1199 : break;
322 0 : case 3:
323 0 : os << ", Zero iterations performed.";
324 0 : break;
325 6 : case 4:
326 6 : os << ", Possible divergence. Peak residual increased by 10% from minimum.";
327 6 : break;
328 7 : case 5:
329 7 : os << ", Exited " << itsAlgorithmName << " minor cycle without reaching any stopping criterion.";
330 7 : break;
331 28 : case 6:
332 28 : os << ", Reached n-sigma threshold.";
333 28 : default:
334 28 : break;
335 : }
336 :
337 3151 : os << LogIO::POST;
338 :
339 3151 : Int rank(0);
340 3151 : String rankStr = EnvironmentVariable::get("OMPI_COMM_WORLD_RANK");
341 3151 : if (!rankStr.empty()) {
342 0 : rank = stoi(rankStr);
343 : }
344 :
345 3151 : int chunkId = chanid; // temporary CAS-13683 workaround
346 : //if (SIMinorCycleController::useSmallSummaryminor()) { // temporary CAS-13683 workaround
347 3151 : if (!fullsummary) { // temporary CAS-13683 workaround
348 3105 : chunkId = chanid + nSubChans*polid;
349 : }
350 3151 : loopcontrols.addSummaryMinor( deconvolverid, chunkId, polid, cycleStartIteration,
351 : startiteration, startmodelflux, startpeakresidual, startpeakresidualnomask,
352 : modelflux, peakresidual, peakresidualnomask, masksum, rank, stopCode, fullsummary);
353 :
354 3151 : loopcontrols.resetCycleIter();
355 :
356 3151 : if( peakresidual > maxResidualAcrossPlanes && stopCode!=0 )
357 1285 : {maxResidualAcrossPlanes=peakresidual; maxResChan=chanid; maxResPol=polid;}
358 :
359 3151 : totalFluxAcrossPlanes += modelflux;
360 :
361 3165 : }// end of polid loop
362 :
363 : }// end of chanid loop
364 :
365 926 : loopcontrols.setPeakResidual( maxResidualAcrossPlanes );
366 :
367 : /// Print total flux over all planes (and max res over all planes). IFF there are more than one plane !!
368 926 : if( nSubChans>1 || nSubPols>1 )
369 : {
370 272 : os << "[" << imagestore->getName() << "] ";
371 272 : os << "Total model flux (over all planes) : " << totalFluxAcrossPlanes; //<< LogIO::POST;
372 272 : os << " Peak Residual (over all planes) : " << maxResidualAcrossPlanes << " in C"<<maxResChan << ":P"<<maxResPol << LogIO::POST;
373 : }
374 :
375 934 : }// end of deconvolve
376 :
377 5988 : Int SDAlgorithmBase::checkStop( SIMinorCycleController &loopcontrols,
378 : Float currentresidual )
379 : {
380 5988 : return loopcontrols.majorCycleRequired(currentresidual);
381 : }
382 :
383 566 : Long SDAlgorithmBase::estimateRAM(const vector<int>& imsize){
384 566 : Long mem=0;
385 566 : if(itsImages)
386 : //Simple deconvolvers will have psf + residual + model + mask (1 plane at a time)
387 0 : mem=sizeof(Float)*(itsImages->getShape()(0))*(itsImages->getShape()(1))*4/1024;
388 566 : else if(imsize.size() >1)
389 566 : mem=sizeof(Float)*(imsize[0])*(imsize[1])*4/1024;
390 : else
391 0 : return 0;
392 : //throw(AipsError("Deconvolver cannot estimate the memory usage at this point"));
393 566 : return mem;
394 : }
395 :
396 1063 : void SDAlgorithmBase::setRestoringBeam( GaussianBeam restbeam, String usebeam )
397 : {
398 2126 : LogIO os( LogOrigin("SDAlgorithmBase","setRestoringBeam",WHERE) );
399 1063 : itsRestoringBeam = restbeam;
400 1063 : itsUseBeam = usebeam;
401 1063 : }
402 :
403 : /*
404 : void SDAlgorithmBase::setMaskOptions( String maskstring )
405 : {
406 : itsMaskString = maskstring;
407 : }
408 : */
409 : /*
410 : void SDAlgorithmBase::loadMask()
411 : {
412 : if (! itsIsMaskLoaded ) {
413 : if( itsMaskString.length()==0 ) {
414 : itsMaskHandler.resetMask( itsImages );
415 : }
416 : else {
417 : itsMaskHandler.fillMask( itsImages->mask() , itsMaskString );
418 : }
419 : itsIsMaskLoaded=true;
420 : }
421 : }
422 : */
423 :
424 739 : void SDAlgorithmBase::restore(std::shared_ptr<SIImageStore> imagestore )
425 : {
426 :
427 1478 : LogIO os( LogOrigin("SDAlgorithmBase","restore",WHERE) );
428 :
429 739 : os << "[" << imagestore->getName() << "] : Restoring model image." << LogIO::POST;
430 :
431 739 : if( imagestore->hasResidualImage() ) imagestore->restore( itsRestoringBeam, itsUseBeam );
432 4 : else{os << "Cannot restore with a residual image" << LogIO::POST;}
433 :
434 :
435 :
436 739 : }
437 :
438 :
439 43 : void SDAlgorithmBase::pbcor(std::shared_ptr<SIImageStore> imagestore )
440 : {
441 :
442 86 : LogIO os( LogOrigin("SDAlgorithmBase","pbcor",WHERE) );
443 :
444 43 : os << "[" << imagestore->getName() << "] : Applying PB correction" << LogIO::POST;
445 :
446 43 : imagestore->pbcor();
447 :
448 43 : }
449 :
450 : /*
451 : void SDAlgorithmBase::restorePlane()
452 : {
453 :
454 : LogIO os( LogOrigin("SDAlgorithmBase","restorePlane",WHERE) );
455 : // << ". Optionally, PB-correct too." << LogIO::POST;
456 :
457 : try
458 : {
459 : // Fit a Gaussian to the PSF.
460 : GaussianBeam beam = itsImages->getPSFGaussian();
461 :
462 : os << "Restore with beam : " << beam.getMajor(Unit("arcmin")) << " arcmin, " << beam.getMinor(Unit("arcmin"))<< " arcmin, " << beam.getPA(Unit("deg")) << " deg)" << LogIO::POST;
463 :
464 : // Initialize restored image
465 : itsImage.set(0.0);
466 : // Copy model into it
467 : itsImage.copyData( LatticeExpr<Float>(itsModel ) );
468 : // Smooth model by beam
469 : StokesImageUtil::Convolve( itsImage, beam);
470 : // Add residual image
471 : itsImage.copyData( LatticeExpr<Float>( itsImage + itsResidual ) );
472 : }
473 : catch(AipsError &x)
474 : {
475 : throw( AipsError("Restoration Error " + x.getMesg() ) );
476 : }
477 :
478 : }
479 : */
480 :
481 : // Use this decide how to partition
482 : // the image for separate calls to 'deconvolve'.
483 : // Give codes to signal one or more of the following.
484 : /// - channel planes separate
485 : /// - stokes planes separate
486 : /// - partitioned-image clean (facets ?)
487 : // Later, can add more complex partitioning schemes....
488 : // but there will be one place to do it, per deconvolver.
489 1001 : void SDAlgorithmBase::queryDesiredShape(Int &nchanchunks, Int &npolchunks, IPosition imshape) // , nImageFacets.
490 : {
491 1001 : AlwaysAssert( imshape.nelements()==4, AipsError );
492 1001 : nchanchunks=imshape(3); // Each channel should appear separately.
493 1001 : npolchunks=imshape(2); // Each pol should appear separately.
494 1001 : }
495 :
496 : /*
497 : /// Make a list of Slices, to send sequentially to the deconvolver.
498 : /// Loop over this list of reference subimages in the 'deconvolve' call.
499 : /// This will support...
500 : /// - channel cube clean
501 : /// - stokes cube clean
502 : /// - partitioned-image clean (facets ?)
503 : /// - 3D deconvolver
504 : void SDAlgorithmBase::partitionImages( std::shared_ptr<SIImageStore> &imagestore )
505 : {
506 : LogIO os( LogOrigin("SDAlgorithmBase","partitionImages",WHERE) );
507 :
508 : IPosition imshape = imagestore->getShape();
509 :
510 :
511 : // TODO : Check which axes is which first !!!
512 : ///// chanAxis_p=CoordinateUtil::findSpectralAxis(dirty_p->coordinates());
513 : //// Vector<Stokes::StokesTypes> whichPols;
514 : //// polAxis_p=CoordinateUtil::findStokesAxis(whichPols, dirty_p->coordinates());
515 : uInt nx = imshape[0];
516 : uInt ny = imshape[1];
517 : uInt npol = imshape[2];
518 : uInt nchan = imshape[3];
519 :
520 : /// (1) /// Set up the Deconvolver Slicers.
521 :
522 : // Ask the deconvolver what shape it wants.
523 : Bool onechan=false, onepol=false;
524 : queryDesiredShape(onechan, onepol);
525 :
526 : uInt nSubImages = ( (onechan)?imshape[3]:1 ) * ( (onepol)?imshape[2]:1 ) ;
527 : uInt polstep = (onepol)?1:npol;
528 : uInt chanstep = (onechan)?1:nchan;
529 :
530 : itsDecSlices.resize( nSubImages );
531 :
532 : uInt subindex=0;
533 : for(uInt pol=0; pol<npol; pol+=polstep)
534 : {
535 : for(uInt chan=0; chan<nchan; chan+=chanstep)
536 : {
537 : AlwaysAssert( subindex < nSubImages , AipsError );
538 : IPosition substart(4,0,0,pol,chan);
539 : IPosition substop(4,nx-1,ny-1, pol+polstep-1, chan+chanstep-1);
540 : itsDecSlices[subindex] = Slicer(substart, substop, Slicer::endIsLast);
541 : subindex++;
542 : }
543 : }
544 :
545 : }// end of partitionImages
546 : */
547 :
548 : /*
549 : void SDAlgorithmBase::initializeSubImages( std::shared_ptr<SIImageStore> &imagestore, uInt subim)
550 : {
551 : itsResidual = SubImage<Float>( *(imagestore->residual()), itsDecSlices[subim], true );
552 : itsPsf = SubImage<Float>( *(imagestore->psf()), itsDecSlices[subim], true );
553 : itsModel = SubImage<Float>( *(imagestore->model()), itsDecSlices[subim], true );
554 :
555 : itsImages = imagestore;
556 :
557 : }
558 : */
559 :
560 : /////////// Helper Functions for all deconvolvers to use if they need it.
561 :
562 0 : Bool SDAlgorithmBase::findMaxAbs(const Array<Float>& lattice,
563 : Float& maxAbs,
564 : IPosition& posMaxAbs)
565 : {
566 : // cout << "findmax : lat shape : " << lattice.shape() << " posmax : " << posMaxAbs << endl;
567 0 : posMaxAbs = IPosition(lattice.shape().nelements(), 0);
568 0 : maxAbs=0.0;
569 : Float minVal;
570 0 : IPosition posmin(lattice.shape().nelements(), 0);
571 0 : minMax(minVal, maxAbs, posmin, posMaxAbs, lattice);
572 : //cout << "min " << minVal << " " << maxAbs << " " << max(lattice) << endl;
573 0 : if(abs(minVal) > abs(maxAbs)){
574 0 : maxAbs=abs(minVal);
575 0 : posMaxAbs=posmin;
576 : }
577 0 : return true;
578 0 : }
579 :
580 2497 : Bool SDAlgorithmBase::findMaxAbsMask(const Array<Float>& lattice,
581 : const Array<Float>& mask,
582 : Float& maxAbs,
583 : IPosition& posMaxAbs)
584 : {
585 :
586 : //cout << "maxabsmask shapes : " << lattice.shape() << " " << mask.shape() << endl;
587 :
588 2497 : posMaxAbs = IPosition(lattice.shape().nelements(), 0);
589 2497 : maxAbs=0.0;
590 : Float minVal;
591 2497 : IPosition posmin(lattice.shape().nelements(), 0);
592 2497 : minMaxMasked(minVal, maxAbs, posmin, posMaxAbs, lattice,mask);
593 : //cout << "min " << minVal << " " << maxAbs << " " << max(lattice) << endl;
594 2497 : if(abs(minVal) > abs(maxAbs)){
595 5 : maxAbs=abs(minVal);
596 5 : posMaxAbs=posmin;
597 : }
598 2497 : return true;
599 2497 : }
600 :
601 : /*
602 : Bool SDAlgorithmBase::findMinMaxMask(const Array<Float>& lattice,
603 : const Array<Float>& mask,
604 : Float& minVal, Float& maxVal)
605 : // IPosition& minPos, IPosition& maxPos)
606 : {
607 : //posMaxAbs = IPosition(lattice.shape().nelements(), 0);
608 : //maxAbs=0.0;
609 : IPosition posmin(lattice.shape().nelements(), 0);
610 : IPosition posmax(lattice.shape().nelements(), 0);
611 : minMaxMasked(minVal, maxVal, posmin, posmax, lattice,mask);
612 :
613 : return true;
614 : }
615 : */
616 : /*
617 :
618 : GaussianBeam SDAlgorithmBase::getPSFGaussian()
619 : {
620 :
621 : GaussianBeam beam;
622 : try
623 : {
624 : if( itsPsf.ndim() > 0 )
625 : {
626 : StokesImageUtil::FitGaussianPSF( itsPsf, beam );
627 : }
628 : }
629 : catch(AipsError &x)
630 : {
631 : LogIO os( LogOrigin("SDAlgorithmBase","getPSFGaussian",WHERE) );
632 : os << "Error in fitting a Gaussian to the PSF : " << x.getMesg() << LogIO::POST;
633 : throw( AipsError("Error in fitting a Gaussian to the PSF" + x.getMesg()) );
634 : }
635 :
636 : return beam;
637 : }
638 :
639 : */
640 : } //# NAMESPACE CASA - END
641 :
|
