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