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