Line data Source code
1 : //# MFCEMemImageSkyModel.cc: Implementation of MFCEMemImageSkyModel class
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/Arrays/Matrix.h>
30 : #include <synthesis/MeasurementComponents/MFCEMemImageSkyModel.h>
31 : #include <casacore/images/Images/PagedImage.h>
32 : #include <casacore/casa/OS/File.h>
33 : #include <casacore/lattices/LRegions/LCBox.h>
34 : #include <casacore/lattices/Lattices/SubLattice.h>
35 : #include <casacore/lattices/Lattices/LatticeStepper.h>
36 : #include <casacore/lattices/Lattices/LatticeUtilities.h>
37 : #include <casacore/lattices/Lattices/LatticeIterator.h>
38 : #include <casacore/lattices/LEL/LatticeExpr.h>
39 : #include <synthesis/MeasurementEquations/SkyEquation.h>
40 : #include <casacore/casa/Exceptions/Error.h>
41 : #include <casacore/casa/BasicSL/String.h>
42 : #include <casacore/casa/Utilities/Assert.h>
43 :
44 : #include <sstream>
45 :
46 : #include <casacore/casa/Logging/LogMessage.h>
47 : #include <casacore/casa/Logging/LogSink.h>
48 : #include <casacore/casa/Logging/LogIO.h>
49 :
50 : #include <casacore/casa/System/Choice.h>
51 :
52 : #include <msvis/MSVis/StokesVector.h>
53 : #include <synthesis/MeasurementEquations/LatConvEquation.h>
54 : #include <synthesis/MeasurementEquations/IncCEMemModel.h>
55 : #include <synthesis/MeasurementEquations/CEMemProgress.h>
56 :
57 : using namespace casacore;
58 : namespace casa {
59 :
60 : // Constructor
61 :
62 0 : MFCEMemImageSkyModel::
63 : MFCEMemImageSkyModel(Float sigma,
64 : Float targetFlux,
65 : Bool constrainFlux,
66 : const Vector<String>& priors,
67 0 : const String& entropy) :
68 0 : CEMemImageSkyModel(sigma, targetFlux, constrainFlux, priors, entropy)
69 : {
70 :
71 : // Hey, we also need some info which controls the change between major cycles!
72 :
73 0 : }
74 :
75 :
76 :
77 : // Clean solver
78 0 : Bool MFCEMemImageSkyModel::solve(SkyEquation& se) {
79 :
80 0 : LogIO os(LogOrigin("MFCEMemImageSkyModel","solve"));
81 :
82 : //Make the PSFs, one per field
83 0 : makeApproxPSFs(se);
84 :
85 : // Validate PSFs for each field
86 0 : Vector<Float> psfmax(numberOfModels()); psfmax=0.0;
87 0 : Vector<Float> psfmaxouter(numberOfModels()); psfmaxouter=0.0;
88 0 : Vector<Float> psfmin(numberOfModels()); psfmin=1.0;
89 0 : Vector<Float> resmax(numberOfModels());
90 0 : Vector<Float> resmin(numberOfModels());
91 :
92 : PagedImage<Float>* priorImagePtr;
93 0 : priorImagePtr=0;
94 0 : if((itsPrior.nelements()>0)&&(itsPrior(0)!="")) {
95 0 : os << "Using prior " << itsPrior(0) << LogIO::POST;
96 0 : priorImagePtr=new PagedImage<Float>(itsPrior(0));
97 : }
98 :
99 0 : Float maxSidelobe=0.0;
100 : Int model;
101 0 : for (model=0;model<numberOfModels();model++) {
102 0 : if(isSolveable(model)) {
103 :
104 0 : IPosition blc(PSF(model).shape().nelements(), 0);
105 0 : IPosition trc(PSF(model).shape()-1);
106 0 : blc(2) = 0; trc(2) = 0;
107 0 : blc(3) = 0; trc(3) = 0;
108 :
109 :
110 0 : SubLattice<Float> subPSF;
111 0 : Int k =0;
112 0 : Int numchan= PSF(model).shape()(3);
113 : //PSF of first non zero plane
114 0 : while(psfmax(model) < 0.1 && k< numchan){
115 0 : blc(3)= k;
116 0 : trc(3)=k;
117 0 : LCBox onePlane(blc, trc, PSF(model).shape());
118 0 : subPSF=SubLattice<Float> ( PSF(model), onePlane, true);
119 : {
120 0 : LatticeExprNode node = max(subPSF);
121 0 : psfmax(model) = node.getFloat();
122 0 : }
123 0 : ++k;
124 0 : }
125 : {
126 0 : LatticeExprNode node = min(subPSF);
127 0 : psfmin(model) = node.getFloat();
128 0 : }
129 : // 32 pixels: pretty arbitrary, but only look for sidelobes
130 : // outside the inner (2n+1) * (2n+1) square
131 0 : psfmaxouter(model) = maxOuter(subPSF, 32);
132 :
133 : os << "Model " << model+1 << ": max, min, maxOuter PSF = "
134 0 : << psfmax(model) << ", " << psfmin(model) << ", " <<
135 0 : psfmaxouter(model) << endl;
136 0 : if(abs(psfmin(model))>maxSidelobe) maxSidelobe=abs(psfmin(model));
137 0 : if(psfmaxouter(model) > maxSidelobe) maxSidelobe= psfmaxouter(model );
138 0 : }
139 : }
140 0 : os << LogIO::POST;
141 :
142 0 : Float absmax=threshold();
143 0 : Float cycleThreshold=0.0;
144 0 : Block< Matrix<Int> > iterations(numberOfModels());
145 0 : Int maxIterations=0;
146 :
147 :
148 0 : Int mpol=image(0).shape()(2);
149 0 : Int mchan=image(0).shape()(3);
150 :
151 0 : Block< Matrix<Float> > alpha(numberOfModels());
152 0 : Block< Matrix<Float> > beta(numberOfModels());
153 0 : Block< Matrix<Float> > Q(numberOfModels());
154 :
155 : // Loop over major cycles
156 0 : Int cycle=0;
157 0 : Bool stop=false;
158 0 : Bool ask=true;
159 0 : Bool modified=false;
160 :
161 0 : if (displayProgress_p) {
162 0 : itsProgress = new CEMemProgress( pgplotter_p );
163 : }
164 :
165 0 : while(absmax>=threshold()&&maxIterations<numberIterations()&&!stop) {
166 :
167 0 : os << "*** Starting major cycle " << cycle+1 << LogIO::POST;
168 0 : cycle++;
169 :
170 : // Make the residual images. We do an incremental update
171 : // for cycles after the first one. If we have only one
172 : // model then we use convolutions to speed the processing
173 0 : Bool incremental(cycle>1);
174 0 : if (!incremental||(itsSubAlgorithm == "full")) {
175 : os << "Using visibility-subtraction for residual calculation"
176 0 : << LogIO::POST;
177 0 : makeNewtonRaphsonStep(se, false);
178 : }
179 : else {
180 : os << "Using XFR-based shortcut for residual calculation"
181 0 : << LogIO::POST;
182 0 : makeNewtonRaphsonStep(se, true);
183 : }
184 : os << "Finished update of residuals"
185 0 : << LogIO::POST;
186 :
187 0 : absmax=maxField(resmax, resmin);
188 :
189 0 : for (model=0;model<numberOfModels();model++) {
190 : os << "Model " << model+1 << ": max, min residuals = "
191 0 : << resmax(model) << ", " << resmin(model) << endl;
192 : }
193 0 : os << LogIO::POST;
194 :
195 : // Can we stop?
196 0 : if(absmax<threshold()) {
197 0 : os << "Reached stopping peak residual = " << absmax << LogIO::POST;
198 0 : stop=true;
199 : }
200 : else {
201 :
202 : // Calculate the threshold for this cycle. Add a safety factor
203 : // This will be fixed someday by an option for an increasing threshold
204 0 : Float fudge = cycleFactor_p * maxSidelobe;
205 0 : if (fudge > 0.8) fudge = 0.8; // painfully slow!
206 :
207 0 : cycleThreshold=max(threshold(), fudge * absmax);
208 : os << "Maximum residual = " << absmax << ", cleaning down to "
209 0 : << cycleThreshold << LogIO::POST;
210 :
211 0 : for (model=0;model<numberOfModels();model++) {
212 :
213 0 : Int npol=image(model).shape()(2);
214 0 : Int nchan=image(model).shape()(3);
215 :
216 0 : IPosition blcDirty(image(model).shape().nelements(), 0);
217 0 : IPosition trcDirty(image(model).shape()-1);
218 :
219 0 : if(cycle==1) {
220 0 : iterations[model].resize(npol, nchan); iterations[model]=0;
221 0 : alpha[model].resize(mpol, mchan); alpha[model] = 0.0;
222 0 : beta[model].resize(mpol, mchan); beta[model] = 0.0;
223 0 : Q[model].resize(mpol, mchan); Q[model] = 0.0;
224 : }
225 :
226 : // Initialize delta image
227 0 : deltaImage(model).set(0.0);
228 :
229 : // Only process solveable models
230 0 : if(isSolveable(model)) {
231 :
232 0 : if((abs(resmax(model))>cycleThreshold)||
233 0 : (abs(resmin(model))>cycleThreshold)) {
234 :
235 0 : os << "Processing model " << model+1 << LogIO::POST;
236 :
237 : // If mask exists, use it;
238 : // If not, use the fluxScale image to figure out
239 0 : Bool doMask = false;
240 0 : Lattice<Float> *maskPointer = 0;
241 0 : if (hasMask(model) && mask(model).nelements() > 1 ) {
242 0 : doMask = true;
243 0 : blcDirty(2)=0; trcDirty(2)=0;
244 0 : blcDirty(3)=0; trcDirty(3)=0;
245 0 : LCBox onePlane(blcDirty, trcDirty, mask(model).shape());
246 :
247 0 : maskPointer = new SubLattice<Float>( mask(model), onePlane, false);
248 :
249 0 : } else if (doFluxScale(model)) {
250 0 : doMask = true;
251 0 : blcDirty(2)=0; trcDirty(2)=0;
252 0 : blcDirty(3)=0; trcDirty(3)=0;
253 0 : LCBox onePlane(blcDirty, trcDirty, fluxScale(model).shape());
254 0 : maskPointer = new SubLattice<Float> ( fluxScale(model), onePlane,
255 0 : true);
256 0 : maskPointer->copyData( (LatticeExpr<Float>)
257 0 : (iif( (*maskPointer > 0.0), 1.0, 0.0) ));
258 0 : }
259 :
260 : // Now deconvolve each channel
261 0 : for (Int chan=0; chan<nchan; chan++) {
262 0 : if(psfmax(model)>0.0) {
263 0 : if(nchan>1) {
264 0 : os<<"Processing channel "<<chan+1<<" of "<<nchan<<LogIO::POST;
265 : }
266 :
267 0 : blcDirty(3) = chan;
268 0 : trcDirty(3) = chan;
269 0 : for (Int pol=0; pol<npol; pol++) {
270 0 : blcDirty(2) = pol; trcDirty(2) = pol;
271 0 : if(npol>1) {
272 0 : os<<"Processing polarization "<<pol+1<<" of "<<npol<<LogIO::POST;
273 : }
274 :
275 0 : LCBox onePlane(blcDirty, trcDirty, image(model).shape());
276 :
277 0 : SubLattice<Float> subImage( image(model), onePlane, true);
278 0 : SubLattice<Float> subResid( residual(model), onePlane);
279 0 : SubLattice<Float> subPSF( PSF(model), onePlane);
280 0 : SubLattice<Float> subDeltaImage( deltaImage(model), onePlane, true);
281 :
282 : // Now make a convolution equation for this
283 : // residual image and psf and then deconvolve it
284 0 : LatConvEquation eqn(subPSF, subResid);
285 :
286 : // Make entropy
287 0 : IncEntropy * myEnt_p = 0;
288 0 : String entString = entropy();
289 0 : if (entString=="mfemptiness") {
290 0 : if (cycle == 1) {
291 0 : os << "Deconvolving with incremental maximum emptiness algorithm" << LogIO::POST;
292 : }
293 0 : myEnt_p = new IncEntropyEmptiness;
294 0 : } else if( entString=="mfentropy") {
295 0 : if (cycle == 1) {
296 0 : os << "Deconvolving with incremental maximum entropy algorithm" << LogIO::POST;
297 : }
298 0 : myEnt_p = new IncEntropyI;
299 : }
300 : else {
301 0 : os << " Known MEM entropies: entropy | emptiness " << LogIO::POST;
302 0 : os << LogIO::SEVERE << "Unknown MEM entropy: " << entString << LogIO::POST;
303 0 : return false;
304 : }
305 :
306 0 : subDeltaImage.set(0.0);
307 : IncCEMemModel *memer;
308 0 : if(priorImagePtr) {
309 0 : if(doMask) {
310 0 : memer= new IncCEMemModel(*myEnt_p, subImage,
311 : subDeltaImage,
312 : *priorImagePtr, *maskPointer,
313 0 : numberIterations(), itsSigma,
314 0 : abs(itsTargetFlux), false,
315 0 : true, false );
316 : }
317 : else {
318 0 : memer=new IncCEMemModel(*myEnt_p, subImage,
319 : subDeltaImage,
320 : *priorImagePtr,
321 0 : numberIterations(), itsSigma,
322 0 : abs(itsTargetFlux), false,
323 0 : true, false );
324 : }
325 : }
326 : else {
327 0 : if(doMask) {
328 0 : memer= new IncCEMemModel(*myEnt_p, subImage,
329 : subDeltaImage,
330 0 : numberIterations(),
331 : *maskPointer,
332 : itsSigma,
333 0 : abs(itsTargetFlux), false,
334 0 : true, false );
335 : }
336 : else {
337 0 : memer=new IncCEMemModel(*myEnt_p, subImage,
338 : subDeltaImage,
339 0 : numberIterations(), itsSigma,
340 0 : abs(itsTargetFlux), false,
341 0 : true, false );
342 : }
343 : }
344 0 : if (displayProgress_p) {
345 0 : memer->setProgress(*itsProgress);
346 : }
347 0 : if (alpha[model](pol, chan) != 0.0) {
348 0 : memer->setAlpha(alpha[model](pol, chan));
349 0 : memer->setBeta(beta[model](pol, chan));
350 0 : memer->setQ(Q[model](pol, chan));
351 : }
352 0 : memer->setGain(Iterate::gain());
353 0 : memer->setThreshold( cycleThreshold );
354 0 : memer->setThresholdSpeedup( cycleSpeedup_p );
355 0 : memer->setInitialNumberIterations(iterations[model](pol, chan));
356 :
357 0 : memer->solve(eqn);
358 :
359 0 : alpha[model](pol, chan) = memer->getAlpha();
360 0 : beta[model](pol, chan) = memer->getBeta();
361 0 : Q[model](pol, chan) = memer->getQ();
362 :
363 0 : iterations[model](pol, chan)=memer->numberIterations();
364 0 : maxIterations=(iterations[model](pol, chan)>maxIterations) ?
365 0 : iterations[model](pol, chan) : maxIterations;
366 0 : modified=true;
367 :
368 0 : if(memer) delete memer; memer=0;
369 :
370 0 : subImage.copyData((LatticeExpr<Float>)
371 0 : (subImage + subDeltaImage ) );
372 :
373 0 : os << "MEM used " << iterations[model](pol, chan) << " iterations"
374 : << " to approach a threshold of " << cycleThreshold
375 0 : << LogIO::POST;
376 0 : }
377 : }
378 : }
379 0 : if (maskPointer) delete maskPointer;
380 : } else {
381 : os<<"Skipping model "<<model+1<<" :peak residual below threshold"
382 0 : <<LogIO::POST;
383 : }
384 : }
385 0 : }
386 : // For now if it has converged withing 20% of maxiter we'll stop
387 :
388 0 : if((maxIterations+Int(numberIterations()/20)) >=numberIterations())
389 0 : stop=true;
390 : //===
391 0 : if(maxIterations<numberIterations()&&ask) {
392 0 : Vector<String> choices(3);
393 0 : choices(0)="Continue";
394 0 : choices(1)="Stop Now";
395 0 : choices(2)="Don't ask again";
396 : String choice =
397 : Choice::choice("MEM iteration: do you want to continue or stop?",
398 0 : choices);
399 0 : if (choice==choices(1)) {
400 0 : os << "Multi-field MEM stopped at user request" << LogIO::POST;
401 0 : stop=true;
402 : }
403 0 : else if (choice==choices(2)) {
404 0 : os << "Continuing: won't ask again" << LogIO::POST;
405 0 : ask=false;
406 : }
407 :
408 0 : }
409 0 : if(!modified) {
410 0 : os << "Nothing happened: stopping" << LogIO::POST;
411 0 : stop=true;
412 : }
413 : }
414 : }
415 0 : if (itsProgress) {
416 0 : delete itsProgress;
417 : }
418 :
419 :
420 0 : if(modified) {
421 0 : os << "Finalizing residual images for all fields" << LogIO::POST;
422 0 : makeNewtonRaphsonStep(se, false);
423 0 : Float finalabsmax=maxField(resmax, resmin);
424 :
425 0 : os << "Final maximum residual = " << finalabsmax << LogIO::POST;
426 :
427 0 : for (model=0;model<numberOfModels();model++) {
428 : os << "Model " << model+1 << ": max, min residuals = "
429 0 : << resmax(model) << ", " << resmin(model) << endl;
430 : }
431 : }
432 : else {
433 0 : os << "Residual images for all fields are up-to-date" << LogIO::POST;
434 : }
435 :
436 0 : os << LogIO::POST;
437 :
438 0 : return(true);
439 0 : };
440 :
441 :
442 : // Find maximum residual
443 0 : Float MFCEMemImageSkyModel::maxField(Vector<Float>& imagemax,
444 : Vector<Float>& imagemin) {
445 :
446 0 : LogIO os(LogOrigin("ImageSkyModel","maxField"));
447 :
448 0 : Float absmax=0.0;
449 0 : imagemax=-1e20;
450 0 : imagemin=1e20;
451 :
452 : // Loop over all models
453 0 : for (Int model=0;model<numberOfModels();model++) {
454 :
455 : // Remember that the residual image can be either as specified
456 : // or created specially.
457 0 : ImageInterface<Float>* imagePtr=0;
458 0 : if(residual_p[model]) {
459 0 : imagePtr=residual_p[model];
460 : }
461 : else {
462 0 : imagePtr=(ImageInterface<Float> *)residualImage_p[model];
463 : }
464 0 : AlwaysAssert(imagePtr, AipsError);
465 0 : AlwaysAssert(imagePtr->shape().nelements()==4, AipsError);
466 0 : Int nx=imagePtr->shape()(0);
467 0 : Int ny=imagePtr->shape()(1);
468 0 : Int npol=imagePtr->shape()(2);
469 :
470 0 : AlwaysAssert((npol==1)||(npol==2)||(npol==4), AipsError);
471 :
472 : // Loop over all channels
473 0 : IPosition onePlane(4, nx, ny, 1, 1);
474 0 : LatticeStepper ls(imagePtr->shape(), onePlane, IPosition(4, 0, 1, 2, 3));
475 0 : RO_LatticeIterator<Float> imageli(*imagePtr, ls);
476 :
477 : // If we are using a mask then reset the region to be
478 : // deconvolved
479 0 : Array<Float> maskArray;
480 :
481 0 : if(hasMask(model)) {
482 0 : LatticeStepper mls(mask(model).shape(), onePlane,
483 0 : IPosition(4, 0, 1, 2, 3));
484 :
485 0 : RO_LatticeIterator<Float> maskli(mask(model), mls);
486 0 : maskli.reset();
487 0 : if (maskli.cursor().shape().nelements() > 1) maskArray=maskli.cursor();
488 0 : }
489 :
490 0 : Int chan=0;
491 : Float imax, imin;
492 0 : imax=-1E20; imagemax(model)=imax;
493 0 : imin=+1E20; imagemin(model)=imin;
494 0 : imageli.reset();
495 :
496 0 : for (imageli.reset();!imageli.atEnd();imageli++,chan++) {
497 :
498 0 : IPosition imageposmax(imageli.cursor().ndim());
499 0 : IPosition imageposmin(imageli.cursor().ndim());
500 :
501 : // If we are using a mask then multiply by it
502 0 : if (hasMask(model)) {
503 0 : Array<Float> limage=imageli.cursor();
504 0 : limage*=maskArray;
505 0 : minMax(imin, imax, imageposmin, imageposmax, limage);
506 0 : }
507 : else {
508 0 : minMax(imin, imax, imageposmin, imageposmax, imageli.cursor());
509 : }
510 0 : if(abs(imax)>absmax) absmax=abs(imax);
511 0 : if(abs(imin)>absmax) absmax=abs(imin);
512 0 : if(imin<imagemin(model)) imagemin(model)=imin;
513 0 : if(imax>imagemax(model)) imagemax(model)=imax;
514 0 : }
515 0 : }
516 0 : return absmax;
517 0 : };
518 :
519 :
520 :
521 : // Find maximum residual
522 0 : Float MFCEMemImageSkyModel::maxOuter(Lattice<Float> & lat, const uInt nCenter )
523 : {
524 0 : Float myMax=0.0;
525 0 : Float myMin=0.0;
526 :
527 : /*
528 : TempLattice<Float> mask(lat.shape());
529 : mask.set(1.0);
530 :
531 : IPosition pos(4,0,0,0,0 );
532 : uInt nxc = lat.shape()(0)/2;
533 : uInt nyc = lat.shape()(1)/2;
534 : for (uInt ix = -nCenter; ix < nCenter; ix++) {
535 : for (uInt iy = -nCenter; iy < nCenter; iy++) {
536 : pos(0) = nxc + ix;
537 : pos(1) = nyc + iy;
538 : mask.putAt( 0.0f, pos ); // mask out the inner section
539 : }
540 : }
541 : {
542 : LatticeExpr<Float> exp = (lat * mask);
543 : LatticeExprNode LEN = max( exp );
544 : myMax = LEN.getFloat();
545 : }
546 : {
547 : LatticeExpr<Float> exp = (lat * mask);
548 : LatticeExprNode LEN = min( exp );
549 : myMin = LEN.getFloat();
550 : }
551 :
552 : */
553 :
554 0 : Array<Float> arr = lat.get();
555 0 : IPosition pos( arr.shape() );
556 0 : uInt nx = lat.shape()(0);
557 0 : uInt ny = lat.shape()(1);
558 0 : uInt nxc = 0;
559 0 : uInt nyc = 0;
560 0 : Float amax = 0.0;
561 0 : for (uInt ix = 0; ix < nx; ix++) {
562 0 : for (uInt iy = 0; iy < ny; iy++) {
563 0 : if (arr(IPosition(4, ix, iy, 0, 0)) > amax) {
564 0 : nxc = ix;
565 0 : nyc = iy;
566 0 : amax = arr(IPosition(4, ix, iy, 0, 0));
567 : }
568 : }
569 : }
570 :
571 0 : uInt nxL = nxc - nCenter;
572 0 : uInt nxH = nxc + nCenter;
573 0 : uInt nyL = nyc - nCenter;
574 0 : uInt nyH = nyc + nCenter;
575 :
576 0 : for (uInt ix = 0; ix < nx; ix++) {
577 0 : for (uInt iy = 0; iy < ny; iy++) {
578 0 : if ( !(ix >= nxL && ix <= nxH && iy >= nyL && iy <= nyH) ) {
579 0 : if (arr(IPosition(4, ix, iy, 0, 0)) > myMax)
580 0 : myMax = arr(IPosition(4, ix, iy, 0, 0));
581 0 : if (arr(IPosition(4, ix, iy, 0, 0)) < myMin)
582 0 : myMin = arr(IPosition(4, ix, iy, 0, 0));
583 : }
584 : }
585 : }
586 :
587 0 : Float absMax = max( abs(myMin), myMax );
588 0 : return absMax;
589 0 : };
590 :
591 : } //#End casa namespace
|
