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: 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 : #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 0 : SDAlgorithmBase::SDAlgorithmBase():
64 0 : itsAlgorithmName("Test")
65 : // itsDecSlices (),
66 : // itsResidual(), itsPsf(), itsModel()
67 : {
68 0 : }
69 :
70 0 : SDAlgorithmBase::~SDAlgorithmBase()
71 : {
72 :
73 0 : }
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 = peakresidualnomask;
140 0 : modelflux = itsImages->getModelFlux();
141 :
142 0 : startpeakresidual = peakresidual;
143 0 : startpeakresidualnomask = peakresidualnomask;
144 0 : startmodelflux = modelflux;
145 :
146 : //Float nsigma = 150.0; // will set by user, fixed for 3sigma for now.
147 0 : Float nsigma = loopcontrols.getNsigma();
148 : //os<<"robustrms nelements="<<robustrms.nelements()<<LogIO::POST;
149 : Float nsigmathresh;
150 0 : if (robustrms.nelements()==0) {
151 0 : nsigmathresh = 0.0;
152 : } else{
153 0 : nsigmathresh = nsigma * (Float)robustrms(IPosition(1,0));
154 : }
155 :
156 : Float thresholdtouse;
157 0 : 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 0 : Array<Double> medians;
161 : //os<<"robuststats rec size="<<robuststats.nfields()<<LogIO::POST;
162 0 : Bool statsexists = false;
163 0 : IPosition statsindex;
164 0 : if (robuststats.nfields()) {
165 : // use existing stats
166 0 : if (robuststats.isDefined("robustrms")) {
167 0 : robuststats.get(RecordFieldId("robustrms"), robustrms);
168 0 : robuststats.get(RecordFieldId("median"), medians);
169 0 : 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 : }
180 : //statsindex=IPosition(1,chanid); // this only support for npol =1, need to fix this
181 0 : if((robustrms.shape().nelements() ==2) && (robustrms.shape()[0] > polid) && (robustrms.shape()[1] > chanid) )
182 0 : statsindex=IPosition(2,polid,chanid);
183 0 : else if((robustrms.shape().nelements() ==1) && (robustrms.shape()[0] > chanid))
184 0 : statsindex=IPosition(1,chanid);
185 : else ///no idea what got passed in the record
186 0 : statsexists=False;
187 :
188 : }
189 0 : if (statsexists) {
190 0 : 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 0 : if (isautomasking) { // new threshold defination
197 : //nsigmathresh = (Float)medians(IPosition(1,0)) + nsigma * (Float)robustrms(IPosition(1,0));
198 0 : nsigmathresh = (Float)medians(statsindex) + nsigma * (Float)robustrms(statsindex);
199 : }
200 : else {
201 0 : nsigmathresh = nsigma * (Float)robustrms(statsindex);
202 : }
203 0 : thresholdtouse = max( nsigmathresh, loopcontrols.getCycleThreshold());
204 : }
205 : else {
206 0 : thresholdtouse = loopcontrols.getCycleThreshold();
207 : }
208 : //os << LogIO::DEBUG1<<"loopcontrols.getCycleThreshold()="<<loopcontrols.getCycleThreshold()<<LogIO::POST;
209 0 : os << LogIO::NORMAL3<<"current CycleThreshold="<<loopcontrols.getCycleThreshold()<<" nsigma threshold="<<nsigmathresh<<LogIO::POST;
210 0 : String thresholddesc = (thresholdtouse == loopcontrols.getCycleThreshold() ? "cyclethreshold" : "n-sigma");
211 0 : os << LogIO::NORMAL3<< "thresholdtouse="<< thresholdtouse << "("<<thresholddesc<<")"<< LogIO::POST;
212 :
213 0 : if (thresholddesc=="n-sigma") {
214 : //os << LogIO::DEBUG1<< "Set nsigma thresh="<<nsigmathresh<<LogIO::POST;
215 0 : loopcontrols.setNsigmaThreshold(nsigmathresh);
216 : }
217 0 : loopcontrols.setPeakResidual( peakresidual );
218 0 : loopcontrols.resetMinResidual(); // Set it to current initial peakresidual.
219 :
220 0 : stopCode = checkStop( loopcontrols, peakresidual );
221 :
222 0 : if( validMask && stopCode==0 )
223 : {
224 : // Init the deconvolver
225 : //where we write in model and residual may be
226 : {
227 :
228 0 : LatticeLocker lockresid (*(itsImages->residual()), FileLocker::Read);
229 0 : LatticeLocker lockmodel (*(itsImages->model()), FileLocker::Read);
230 0 : LatticeLocker lockmask (*(itsImages->mask()), FileLocker::Read);
231 0 : LatticeLocker lockpsf (*(itsImages->psf()), FileLocker::Read);
232 :
233 0 : initializeDeconvolver();
234 : }
235 :
236 0 : while ( stopCode==0 )
237 : {
238 :
239 0 : if (nsigma>0.0) {
240 0 : os << "Using " << thresholddesc << " for threshold criterion: (cyclethreshold="<<loopcontrols.getCycleThreshold()<< ", nsigma threshold="<<nsigmathresh<<" )" << LogIO::POST;
241 0 : loopcontrols.setNsigmaThreshold(nsigmathresh);
242 : }
243 0 : Int thisniter = loopcontrols.getCycleNiter() <5000 ? loopcontrols.getCycleNiter() : 2000;
244 :
245 0 : 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 0 : takeOneStep( loopcontrols.getLoopGain(),
254 : // loopcontrols.getCycleNiter(),
255 : thisniter,
256 : //loopcontrols.getCycleThreshold(),
257 : thresholdtouse,
258 : peakresidual,
259 : modelflux,
260 0 : iterdone);
261 : }
262 :
263 0 : os << LogIO::NORMAL1 << "SDAlgoBase: After one step, dec : " << deconvolverid << " residual=" << peakresidual << " model=" << modelflux << " iters=" << iterdone << LogIO::POST;
264 :
265 0 : SynthesisUtilMethods::getResource("In Deconvolver : one step" );
266 :
267 0 : loopcontrols.incrementMinorCycleCount( iterdone ); // CAS-8767 : add subimageindex and merge with addSummaryMinor call later.
268 :
269 0 : 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 0 : if(stopCode==0 && iterdone != thisniter)
277 : {
278 0 : os << LogIO::NORMAL1 << "Exited " << itsAlgorithmName << " minor cycle without satisfying stopping criteria " << LogIO::POST;
279 0 : 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 0 : LatticeLocker lock1 (*(itsImages->residual()), FileLocker::Write);
288 0 : LatticeLocker lock2 (*(itsImages->model()), FileLocker::Write);
289 0 : finalizeDeconvolver();
290 : }
291 :
292 : // get the new peakres without a mask for the summary minor
293 0 : peakresidualnomask = itsImages->getPeakResidual();
294 :
295 : }// if validmask
296 :
297 : // same as checking on itscycleniter.....
298 0 : loopcontrols.setUpdatedModelFlag( loopcontrols.getIterDone()-startiteration );
299 :
300 0 : os << "[" << imagestore->getName();
301 0 : if(nSubChans>1) os << ":C" << chanid ;
302 0 : if(nSubPols>1) os << ":P" << polid ;
303 0 : 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 0 : << ", peakres=" << startpeakresidual << "->" << peakresidual ;
309 :
310 0 : switch (stopCode)
311 : {
312 0 : case 0:
313 0 : os << ", Skipped this plane. Zero mask.";
314 0 : break;
315 0 : case 1:
316 0 : os << ", Reached cycleniter.";
317 0 : break;
318 0 : case 2:
319 0 : os << ", Reached cyclethreshold.";
320 0 : break;
321 0 : case 3:
322 0 : os << ", Zero iterations performed.";
323 0 : break;
324 0 : case 4:
325 0 : os << ", Possible divergence. Peak residual increased by 10% from minimum.";
326 0 : break;
327 0 : case 5:
328 0 : os << ", Exited " << itsAlgorithmName << " minor cycle without reaching any stopping criterion.";
329 0 : break;
330 0 : case 6:
331 0 : os << ", Reached n-sigma threshold.";
332 0 : default:
333 0 : break;
334 : }
335 :
336 0 : os << LogIO::POST;
337 :
338 0 : Int rank(0);
339 0 : String rankStr = EnvironmentVariable::get("OMPI_COMM_WORLD_RANK");
340 0 : if (!rankStr.empty()) {
341 0 : rank = stoi(rankStr);
342 : }
343 :
344 0 : int chunkId = chanid; // temporary CAS-13683 workaround
345 : //if (SIMinorCycleController::useSmallSummaryminor()) { // temporary CAS-13683 workaround
346 0 : if (!fullsummary) { // temporary CAS-13683 workaround
347 0 : chunkId = chanid + nSubChans*polid;
348 : }
349 0 : loopcontrols.addSummaryMinor( deconvolverid, chunkId, polid, cycleStartIteration,
350 : startiteration, startmodelflux, startpeakresidual, startpeakresidualnomask,
351 : modelflux, peakresidual, peakresidualnomask, masksum, rank, stopCode, fullsummary);
352 :
353 0 : loopcontrols.resetCycleIter();
354 :
355 0 : if( peakresidual > maxResidualAcrossPlanes && stopCode!=0 )
356 0 : {maxResidualAcrossPlanes=peakresidual; maxResChan=chanid; maxResPol=polid;}
357 :
358 0 : totalFluxAcrossPlanes += modelflux;
359 :
360 : }// end of polid loop
361 :
362 : }// end of chanid loop
363 :
364 0 : 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 0 : if( nSubChans>1 || nSubPols>1 )
368 : {
369 0 : os << "[" << imagestore->getName() << "] ";
370 0 : os << "Total model flux (over all planes) : " << totalFluxAcrossPlanes; //<< LogIO::POST;
371 0 : os << " Peak Residual (over all planes) : " << maxResidualAcrossPlanes << " in C"<<maxResChan << ":P"<<maxResPol << LogIO::POST;
372 : }
373 :
374 0 : }// end of deconvolve
375 :
376 0 : Int SDAlgorithmBase::checkStop( SIMinorCycleController &loopcontrols,
377 : Float currentresidual )
378 : {
379 0 : return loopcontrols.majorCycleRequired(currentresidual);
380 : }
381 :
382 0 : Long SDAlgorithmBase::estimateRAM(const vector<int>& imsize){
383 0 : Long mem=0;
384 0 : 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 0 : else if(imsize.size() >1)
388 0 : 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 0 : return mem;
393 : }
394 :
395 0 : void SDAlgorithmBase::setRestoringBeam( GaussianBeam restbeam, String usebeam )
396 : {
397 0 : LogIO os( LogOrigin("SDAlgorithmBase","setRestoringBeam",WHERE) );
398 0 : itsRestoringBeam = restbeam;
399 0 : itsUseBeam = usebeam;
400 0 : }
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 0 : void SDAlgorithmBase::restore(std::shared_ptr<SIImageStore> imagestore )
424 : {
425 :
426 0 : LogIO os( LogOrigin("SDAlgorithmBase","restore",WHERE) );
427 :
428 0 : os << "[" << imagestore->getName() << "] : Restoring model image." << LogIO::POST;
429 :
430 0 : if( imagestore->hasResidualImage() ) imagestore->restore( itsRestoringBeam, itsUseBeam );
431 0 : else{os << "Cannot restore with a residual image" << LogIO::POST;}
432 :
433 :
434 :
435 0 : }
436 :
437 :
438 0 : void SDAlgorithmBase::pbcor(std::shared_ptr<SIImageStore> imagestore )
439 : {
440 :
441 0 : LogIO os( LogOrigin("SDAlgorithmBase","pbcor",WHERE) );
442 :
443 0 : os << "[" << imagestore->getName() << "] : Applying PB correction" << LogIO::POST;
444 :
445 0 : imagestore->pbcor();
446 :
447 0 : }
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 0 : void SDAlgorithmBase::queryDesiredShape(Int &nchanchunks, Int &npolchunks, IPosition imshape) // , nImageFacets.
489 : {
490 0 : AlwaysAssert( imshape.nelements()==4, AipsError );
491 0 : nchanchunks=imshape(3); // Each channel should appear separately.
492 0 : npolchunks=imshape(2); // Each pol should appear separately.
493 0 : }
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 : }
578 :
579 0 : 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 0 : posMaxAbs = IPosition(lattice.shape().nelements(), 0);
588 0 : maxAbs=0.0;
589 : Float minVal;
590 0 : IPosition posmin(lattice.shape().nelements(), 0);
591 0 : minMaxMasked(minVal, maxAbs, posmin, posMaxAbs, lattice,mask);
592 : //cout << "min " << minVal << " " << maxAbs << " " << max(lattice) << endl;
593 0 : if(abs(minVal) > abs(maxAbs)){
594 0 : maxAbs=abs(minVal);
595 0 : posMaxAbs=posmin;
596 : }
597 0 : return true;
598 : }
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 :
|
