Line data Source code
1 : //# Copyright (C) 1998,1999,2000,2001,2003
2 : //# Associated Universities, Inc. Washington DC, USA.
3 : //#
4 : //# This program is free software; you can redistribute it and/or modify it
5 : //# under the terms of the GNU General Public License as published by the Free
6 : //# Software Foundation; either version 2 of the License, or (at your option)
7 : //# any later version.
8 : //#
9 : //# This program is distributed in the hope that it will be useful, but WITHOUT
10 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
12 : //# more details.
13 : //#
14 : //# You should have received a copy of the GNU General Public License along
15 : //# with this program; if not, write to the Free Software Foundation, Inc.,
16 : //# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
17 : //#
18 : //# Correspondence concerning AIPS++ should be addressed as follows:
19 : //# Internet email: aips2-request@nrao.edu.
20 : //# Postal address: AIPS++ Project Office
21 : //# National Radio Astronomy Observatory
22 : //# 520 Edgemont Road
23 : //# Charlottesville, VA 22903-2475 USA
24 : //#
25 :
26 : #ifndef IMAGEANALYSIS_IMAGESTATSCALCULATOR_TCC
27 : #define IMAGEANALYSIS_IMAGESTATSCALCULATOR_TCC
28 :
29 : #include <imageanalysis/ImageAnalysis/ImageStatsCalculator.h>
30 :
31 : #include <casacore/casa/BasicSL/String.h>
32 : #include <casacore/images/Images/ImageUtilities.h>
33 : #include <imageanalysis/ImageAnalysis/SubImageFactory.h>
34 :
35 : #include <iomanip>
36 :
37 : using namespace casacore;
38 : namespace casa {
39 :
40 : template <class T>
41 : const String ImageStatsCalculator<T>::_class = "ImageStatsCalculator";
42 :
43 : template <class T>
44 : const String ImageStatsCalculator<T>::SIGMA = "sigma";
45 :
46 0 : template <class T> ImageStatsCalculator<T>::ImageStatsCalculator(
47 : const SPCIIT image,
48 : const Record *const ®ionPtr,
49 : const String& maskInp,
50 : Bool beVerboseDuringConstruction
51 : ) : ImageStatsBase<T>(
52 : image, regionPtr, maskInp
53 0 : ) {
54 0 : this->_construct(beVerboseDuringConstruction);
55 0 : }
56 :
57 0 : template <class T> ImageStatsCalculator<T>::~ImageStatsCalculator() {}
58 :
59 : template <class T> Record ImageStatsCalculator<T>::calculate() {
60 : *this->_getLog() << LogOrigin(_class, __func__);
61 : std::unique_ptr<std::vector<String> > messageStore(
62 : this->_getLogFile() ? new std::vector<String>() : nullptr
63 : );
64 : Record retval = statistics(messageStore.get());
65 : Bool writeFile = this->_openLogfile();
66 : if (_verbose || writeFile) {
67 : if (writeFile) {
68 : for (
69 : auto iter = messageStore->begin();
70 : iter != messageStore->end(); ++iter
71 : ) {
72 : this->_writeLogfile("# " + *iter, false, false);
73 : }
74 : }
75 : IPosition shape = _axes.empty() ? IPosition(_subImage->ndim(), 1)
76 : : _subImage->shape();
77 : for (const auto& axis: _axes) {
78 : shape[axis] = 1;
79 : }
80 : Record r;
81 : auto csys = _subImage->coordinates();
82 : csys.save(r, "");
83 : try {
84 : auto tempIm = ImageFactory::fromShape<T>(casacore::String(""), shape.asVector(), r);
85 : _reportDetailedStats(tempIm, retval);
86 : }
87 : catch (const AipsError& x) {
88 : *this->_getLog() << LogIO::WARN << "Unable to collapse image "
89 : << "so detailed per plane statistics reporting is not "
90 : << "possible. The exception message was " << x.getMesg()
91 : << LogIO::POST;
92 : }
93 : }
94 : _sanitizeDueToRegionSelection(retval);
95 : return retval;
96 : }
97 :
98 : template <class T> void
99 : ImageStatsCalculator<T>::_sanitizeDueToRegionSelection(Record& retval) const {
100 : if (_axes.empty()) {
101 : return;
102 : }
103 : if (! this->_getRegion() || this->_getRegion()->empty()) {
104 : // no region selection, nothing to sanitize
105 : return;
106 : }
107 : // create subimage template based on region only
108 : TempImage<T> tempIm(this->_getImage()->shape(), this->_getImage()->coordinates());
109 : auto subim = SubImageFactory<T>::createSubImageRO(
110 : tempIm, *this->_getRegion(), "", nullptr, AxesSpecifier(), False
111 : );
112 : if (! subim->isMasked()) {
113 : // no pixels masked because of region selection
114 : return;
115 : }
116 : auto ndim = subim->ndim();
117 : auto allAxes = IPosition::makeAxisPath(ndim);
118 : IPosition cursor;
119 : for (auto a: _axes) {
120 : cursor.append(IPosition(1, a));
121 : }
122 : auto displayAxes = allAxes.otherAxes(ndim, cursor);
123 : // key is axis number, value is set of planes that are completely masked
124 : std::map<uInt, std::set<uInt>> excludePlanes;
125 : Bool mustExclude = False;
126 : for (auto d: displayAxes) {
127 : excludePlanes[d] = std::set<uInt>();
128 : IPosition cursorShape = subim->shape();
129 : cursorShape[d] = 1;
130 : RO_MaskedLatticeIterator<T> lattIter(*subim, cursorShape);
131 : uInt planeNum = 0;
132 : for (lattIter.atStart(); ! lattIter.atEnd(); ++lattIter, ++planeNum) {
133 : if (! anyTrue(lattIter.getMask())) {
134 : excludePlanes[d].insert(planeNum);
135 : mustExclude = True;
136 : }
137 : }
138 : }
139 : if (! mustExclude) {
140 : // no planes to exclude
141 : return;
142 : }
143 : auto nfields = retval.nfields();
144 : // n is the index of the axis within the displayAxes
145 : uInt n = 0;
146 : for (auto d: displayAxes) {
147 : if (excludePlanes[d].empty()) {
148 : // no planes to exclude for this axis
149 : continue;
150 : }
151 : for (uInt i=0; i<nfields; ++i) {
152 : auto fieldName = retval.name(i);
153 : if (fieldName == "blc" || fieldName == "trc") {
154 : continue;
155 : }
156 : if (isArray(retval.dataType(i))) {
157 : switch (retval.dataType(i)) {
158 : case TpArrayDouble: {
159 : auto x = retval.asArrayDouble(i);
160 : _removePlanes(x, n, excludePlanes[d]);
161 : retval.define(i, x);
162 : break;
163 : }
164 : case TpArrayInt: {
165 : auto x = retval.asArrayInt(i);
166 : _removePlanes(x, n, excludePlanes[d]);
167 : retval.define(i, x);
168 : break;
169 : }
170 : default:
171 : ThrowCc("Unhandled data type");
172 : break;
173 : }
174 : }
175 : }
176 : ++n;
177 : }
178 : }
179 :
180 : template <class T> template <class U>
181 : void ImageStatsCalculator<T>::_removePlanes(
182 : Array<U>& arr, uInt axis, const std::set<uInt>& planes
183 : ) const {
184 : IPosition oldShape = arr.shape();
185 : IPosition newShape = oldShape;
186 : newShape[axis] -= planes.size();
187 : Array<U> newArray(newShape);
188 : // do a plane by plane copy into the new array
189 : auto nOldPlanes = oldShape[axis];
190 : auto begin = planes.begin();
191 : auto end = planes.end();
192 : auto ndim = arr.ndim();
193 : IPosition newSliceStart(ndim, 0);
194 : IPosition newSliceEnd = newShape - 1;
195 : newSliceEnd[axis] = 0;
196 : IPosition oldSliceStart(ndim, 0);
197 : IPosition oldSliceEnd = oldShape - 1;
198 : oldSliceEnd[axis] = 0;
199 : Slicer newSlice(newSliceStart, newSliceEnd, Slicer::endIsLast);
200 : Slicer oldSlice(oldSliceStart, oldSliceEnd, Slicer::endIsLast);
201 : for (uInt i=0; i<nOldPlanes; ++i, ++oldSliceStart[axis], ++oldSliceEnd[axis]) {
202 : if (std::find(begin, end, i) == end) {
203 : newSlice.setStart(newSliceStart);
204 : newSlice.setEnd(newSliceEnd);
205 : oldSlice.setStart(oldSliceStart);
206 : oldSlice.setEnd(oldSliceEnd);
207 : newArray(newSlice) = arr(oldSlice);
208 : ++newSliceStart[axis];
209 : ++newSliceEnd[axis];
210 : }
211 : }
212 : arr.assign(newArray);
213 : }
214 :
215 : template <class T> void ImageStatsCalculator<T>::setVerbose(Bool v) {
216 : if (_verbose != v) {
217 : this->_resetStats();
218 : }
219 : _verbose = v;
220 : }
221 :
222 : template <class T> void ImageStatsCalculator<T>::setDisk(Bool d) {
223 : if (_disk != d) {
224 : this->_resetStats();
225 : }
226 : _disk = d;
227 : }
228 :
229 : template <class T> void ImageStatsCalculator<T>::_reportDetailedStats(
230 : const SPCIIT tempIm, const Record& retval
231 : ) {
232 : auto nptsArr = retval.asArrayDouble("npts");
233 : if (nptsArr.empty()) {
234 : auto msg = "NO UNMASKED POINTS FOUND, NO STATISTICS WERE COMPUTED";
235 : *this->_getLog() << LogIO::NORMAL << msg << LogIO::POST;
236 : if (this->_getLogFile()) {
237 : this->_writeLogfile(msg, false, false);
238 : }
239 : this->_closeLogfile();
240 : return;
241 : }
242 : const CoordinateSystem& csys = tempIm->coordinates();
243 : auto worldAxes = csys.worldAxisNames();
244 : auto imShape = tempIm->shape();
245 : vector<uInt> colwidth;
246 : Int stokesCol = -1;
247 : Int freqCol = -1;
248 : Int raCol = -1;
249 : Int decCol = -1;
250 : IPosition otherCol;
251 : for (Int i=worldAxes.size()-1; i>=0; i--) {
252 : auto gg = worldAxes[i];
253 : gg.upcase();
254 : if (gg == "RIGHT ASCENSION") {
255 : raCol = i;
256 : }
257 : else if (gg == "DECLINATION") {
258 : decCol = i;
259 : }
260 : else if (gg == "FREQUENCY") {
261 : freqCol = i;
262 : }
263 : else if (gg == "STOKES") {
264 : stokesCol = i;
265 : }
266 : else {
267 : otherCol.append(IPosition(1, i));
268 : }
269 : }
270 : IPosition idx(worldAxes.size(), 0);
271 : uInt myloc = 0;
272 : IPosition reportAxes;
273 : if (stokesCol >= 0) {
274 : idx[myloc] = stokesCol;
275 : if (imShape[stokesCol] > 1) {
276 : reportAxes.prepend(IPosition(1, stokesCol));
277 : }
278 : ++myloc;
279 : }
280 : if (freqCol >= 0) {
281 : idx[myloc] = freqCol;
282 : if (imShape[freqCol] > 1) {
283 : reportAxes.prepend(IPosition(1, freqCol));
284 : }
285 : ++myloc;
286 : }
287 : if (decCol >= 0) {
288 : idx[myloc] = decCol;
289 : if (imShape[decCol] > 1) {
290 : reportAxes.prepend(IPosition(1, decCol));
291 : }
292 : ++myloc;
293 : }
294 : if (raCol >= 0) {
295 : idx[myloc] = raCol;
296 : if (imShape[raCol] > 1) {
297 : reportAxes.prepend(IPosition(1, raCol));
298 : }
299 : myloc++;
300 : }
301 : if (otherCol.size() > 0) {
302 : for (uInt i=0; i<otherCol.nelements(); ++i) {
303 : idx[myloc] = otherCol[i];
304 : ++myloc;
305 : if (imShape[otherCol[i]] > 1) {
306 : reportAxes.append(IPosition(1, otherCol[i]));
307 : }
308 : }
309 : }
310 : Bool doVelocity = csys.hasSpectralAxis()
311 : && csys.spectralCoordinate().restFrequency() > 0;
312 : ostringstream oss;
313 : // CSSC wants "#" in log file but not in logger output, sigh
314 : for (auto ax : reportAxes) {
315 : if (ax == freqCol) {
316 : if (doVelocity) {
317 : oss << "VELOCITY column unit = "
318 : << csys.spectralCoordinate().velocityUnit() << endl;
319 : }
320 : else {
321 : oss << "FREQUENCY column unit = "
322 : << csys.spectralCoordinate().worldAxisUnits()[0] << endl;
323 : }
324 : if (_verbose) {
325 : *this->_getLog() << LogIO::NORMAL << oss.str() << LogIO::POST;
326 : }
327 : if (this->_getLogFile()) {
328 : this->_writeLogfile("#" + oss.str(), false, false);
329 : }
330 : oss.str("");
331 : }
332 : }
333 : auto bUnit = this->_getImage()->units().getName();
334 : const auto alg = this->_getAlgorithm();
335 : const auto doBiweight = alg == StatisticsData::BIWEIGHT;
336 : if (_verbose) {
337 : oss.str("");
338 : if (! doBiweight) {
339 : oss << "Sum column unit = " << bUnit << endl;
340 : }
341 : oss << "Mean column unit = " << bUnit << endl;
342 : oss << "Std_dev column unit = " << bUnit << endl;
343 : oss << "Minimum column unit = " << bUnit << endl;
344 : oss << "Maximum column unit = " << bUnit << endl;
345 : *this->_getLog() << LogIO::NORMAL << oss.str() << LogIO::POST;
346 : oss.str("");
347 : }
348 : if (this->_getLogFile()) {
349 : oss.str("");
350 : if (! doBiweight) {
351 : oss << "#Sum column unit = " << bUnit << endl;
352 : }
353 : oss << "#Mean column unit = " << bUnit << endl;
354 : oss << "#Std_dev column unit = " << bUnit << endl;
355 : oss << "#Minimum column unit = " << bUnit << endl;
356 : oss << "#Maximum column unit = " << bUnit << endl;
357 : this->_writeLogfile(oss.str(), false, false);
358 : oss.str("");
359 : }
360 : for (auto ax : reportAxes) {
361 : String gg = worldAxes[ax];
362 : gg.upcase();
363 : uInt width = gg == "STOKES" ? 6 : gg == "FREQUENCY"? 16: 15;
364 : if (
365 : gg == "FREQUENCY" && doVelocity
366 : ) {
367 : gg = "VELOCITY";
368 : }
369 : colwidth.push_back(width);
370 : oss << setw(width) << gg << " "
371 : << gg << "(Plane)" << " ";
372 : width = gg.size() + 8;
373 : colwidth.push_back(width);
374 : }
375 : Vector<Int> axesMap = reportAxes.asVector();
376 : GenSort<Int>::sort(axesMap);
377 : if (doBiweight) {
378 : oss << "Npts Mean Std_dev Minimum Maximum ";
379 : }
380 : else {
381 : oss << "Npts Sum Mean Rms Std_dev Minimum Maximum ";
382 : }
383 : std::map<String, uInt> chauvIters;
384 : const auto& stats = this->_getImageStats();
385 : if (alg == StatisticsData::CHAUVENETCRITERION) {
386 : chauvIters = stats->getChauvenetNiter();
387 : oss << " N Iter";
388 : }
389 : oss << endl;
390 : if (_verbose) {
391 : *this->_getLog() << LogIO::NORMAL << oss.str() << LogIO::POST;
392 : }
393 : if (this->_getLogFile()) {
394 : this->_writeLogfile("#" + oss.str(), false, false);
395 : }
396 : oss.str("");
397 : for (uInt i=0; i<7; ++i) {
398 : colwidth.push_back(12);
399 : }
400 : if (alg == StatisticsData::CHAUVENETCRITERION) {
401 : colwidth.push_back(6);
402 : }
403 : TileStepper ts(
404 : tempIm->niceCursorShape(),
405 : IPosition(tempIm->ndim(), 1), idx
406 : );
407 : RO_MaskedLatticeIterator<T> inIter(
408 : *tempIm, ts
409 : );
410 : Vector<Double> world;
411 : IPosition arrayIndex(axesMap.nelements(), 0);
412 : IPosition blc = stats->getBlc();
413 : IPosition position(tempIm->ndim());
414 : uInt width = 13;
415 : Vector<Vector<String> > coords(reportAxes.size());
416 : auto i = 0;
417 : for (const auto& axis: reportAxes) {
418 : Vector<Double> indices = indgen(imShape[axis], 0.0, 1.0);
419 : uInt prec = axis == freqCol ? 9 : 5;
420 : if (doVelocity && reportAxes[i] == freqCol) {
421 : const SpectralCoordinate& spc = csys.spectralCoordinate();
422 : Vector<Double> vels;
423 : spc.pixelToVelocity(vels, indices);
424 : vector<String> sv;
425 : for (const auto& v : vels) {
426 : ostringstream oss;
427 : oss << setprecision(prec) << v;
428 : sv.push_back(oss.str());
429 : }
430 : coords[i] = Vector<String>(sv);
431 : }
432 : else {
433 : ImageUtilities::pixToWorld(
434 : coords[i], csys, axis, _axes,
435 : IPosition(imShape.size(),0), imShape-1, indices, prec,
436 : true
437 : );
438 : }
439 : ++i;
440 : }
441 : uInt count = 0;
442 : for (inIter.reset(); ! inIter.atEnd(); ++inIter) {
443 : oss << std::scientific;
444 : uInt colNum = 0;
445 : position = inIter.position();
446 : csys.toWorld(world, position);
447 : if (axesMap.empty()) {
448 : arrayIndex = IPosition(1, 0);
449 : }
450 : else {
451 : auto n = axesMap.nelements();
452 : for (uInt i=0; i<n; ++i) {
453 : arrayIndex[i] = position[axesMap[i]];
454 : }
455 : }
456 : auto npts = nptsArr(arrayIndex);
457 : if (npts == 0) {
458 : // CAS-10183, do not log planes for which there are no good points
459 : continue;
460 : }
461 : for (uInt i=0; i<reportAxes.nelements(); ++i) {
462 : oss << setw(colwidth[colNum]);
463 : oss << coords[i][position[reportAxes[i]]];
464 : ++colNum;
465 : oss << " " << setw(colwidth[colNum])
466 : << (position[reportAxes[i]] + blc[reportAxes[i]]) << " ";
467 : ++colNum;
468 : }
469 : oss << std::setw(width) << npts << " ";
470 : if (alg != StatisticsData::BIWEIGHT) {
471 : oss << std::setw(width) << retval.asArrayDouble("sum")(arrayIndex) << " ";
472 : }
473 : oss << std::setw(width) << retval.asArrayDouble("mean")(arrayIndex) << " ";
474 : if (alg != StatisticsData::BIWEIGHT) {
475 : oss << std::setw(width) << retval.asArrayDouble("rms")(arrayIndex) << " ";
476 : }
477 : oss << std::setw(width) << retval.asArrayDouble(SIGMA)(arrayIndex) << " "
478 : << std::setw(width) << retval.asArrayDouble("min")(arrayIndex) << " "
479 : << std::setw(width) << retval.asArrayDouble("max")(arrayIndex);
480 : if (alg == StatisticsData::CHAUVENETCRITERION) {
481 : ostringstream pos;
482 : pos << position;
483 : oss << std::setw(6) << " " << chauvIters[pos.str()];
484 : ++count;
485 : }
486 : oss << endl;
487 : if (_verbose) {
488 : *this->_getLog() << LogIO::NORMAL << oss.str() << LogIO::POST;
489 : }
490 : // add a space at the beginning of the line to account for the
491 : // "#" in the column header
492 : this->_writeLogfile(" " + oss.str(), false, false);
493 : oss.str("");
494 : }
495 : this->_closeLogfile();
496 : }
497 :
498 0 : template <class T> Record ImageStatsCalculator<T>::statistics(
499 : std::vector<String> *const &messageStore
500 : ) {
501 0 : LogOrigin myOrigin(_class, __func__);
502 0 : *this->_getLog() << myOrigin;
503 0 : CountedPtr<ImageRegion> region, mask;
504 0 : String mtmp = this->_getMask();
505 0 : if (mtmp == "false" || mtmp == "[]") {
506 0 : mtmp = "";
507 : }
508 0 : _subImage = SubImageFactory<T>::createSubImageRO(
509 0 : region, mask, *this->_getImage(), *this->_getRegion(), mtmp,
510 0 : (_verbose ? this->_getLog().get() : 0), AxesSpecifier(),
511 : this->_getStretch()
512 : );
513 0 : *this->_getLog() << myOrigin;
514 : // Find BLC of _subImage in pixels and world coords, and output the
515 : // information to the logger.
516 : // NOTE: ImageStatitics can't do this because it only gets the _subImage
517 : // not a region and the full image.
518 0 : IPosition shape = _subImage->shape();
519 0 : IPosition blc(_subImage->ndim(), 0);
520 0 : IPosition trc(shape - 1);
521 0 : if (region) {
522 0 : LatticeRegion latRegion = region->toLatticeRegion(
523 0 : this->_getImage()->coordinates(), this->_getImage()->shape()
524 : );
525 0 : Slicer sl = latRegion.slicer();
526 0 : blc = sl.start();
527 0 : trc = sl.end();
528 : }
529 : // for precision
530 0 : CoordinateSystem csys = this->_getImage()->coordinates();
531 0 : Int precis = -1;
532 0 : if (csys.hasDirectionCoordinate()) {
533 0 : DirectionCoordinate dirCoord = csys.directionCoordinate();
534 0 : Vector<String> dirUnits = dirCoord.worldAxisUnits();
535 0 : Vector<Double> dirIncs = dirCoord.increment();
536 0 : for (uInt i=0; i< dirUnits.size(); ++i) {
537 0 : Quantity inc(dirIncs[i], dirUnits[i]);
538 0 : inc.convert("s");
539 0 : Int newPrecis = abs(int(floor(log10(inc.getValue()))));
540 0 : precis = (newPrecis > 2 && newPrecis > precis) ? newPrecis : precis;
541 : }
542 : }
543 0 : String blcf, trcf;
544 0 : blcf = CoordinateUtil::formatCoordinate(blc, csys, precis);
545 0 : trcf = CoordinateUtil::formatCoordinate(trc, csys, precis);
546 0 : auto& stats = this->_getImageStats();
547 0 : if (! stats) {
548 0 : this->_resetStats(
549 0 : _verbose
550 0 : ? new ImageStatistics<T> (*_subImage, *this->_getLog(), true, _disk)
551 0 : : new ImageStatistics<T> (*_subImage, true, _disk)
552 : );
553 : }
554 : else {
555 0 : stats->resetError();
556 0 : if (
557 0 : _haveRegionsChanged(
558 : region.get(), mask.get(),
559 : _oldStatsRegion.get(), _oldStatsMask.get()
560 : )
561 : ) {
562 0 : stats->setNewImage(*_subImage);
563 : }
564 : }
565 : // prevent the table of stats we no longer use from being logged
566 0 : stats->setListStats(false);
567 0 : auto myAlg = this->_configureAlgorithm();
568 0 : _logStartup(
569 : messageStore, myAlg, blc, trc, blcf, trcf
570 : );
571 0 : auto doBiweight = this->_getAlgorithm() == StatisticsData::BIWEIGHT;
572 0 : if (_robust && doBiweight) {
573 0 : *this->_getLog() << LogIO::WARN << "The biweight algorithm does not "
574 : << "support the computation of quantile-like statistics. "
575 0 : << "These will not be computed" << LogIO::POST;
576 0 : _robust = False;
577 : }
578 0 : stats->setComputeQuantiles(_robust);
579 0 : if (messageStore) {
580 0 : stats->recordMessages(true);
581 : }
582 0 : stats->setPrecision(precis);
583 0 : stats->setBlc(blc);
584 : // Assign old regions to current regions
585 0 : _oldStatsMask.reset();
586 0 : _oldStatsRegion = region;
587 0 : _oldStatsMask = mask;
588 : // Set cursor axes
589 0 : *this->_getLog() << myOrigin;
590 0 : ThrowIf(! stats->setAxes(_axes), stats->errorMessage());
591 0 : ThrowIf(
592 : !stats->setInExCludeRange(_includepix, _excludepix, false),
593 : stats->errorMessage()
594 : );
595 : // Tell what to list
596 0 : ThrowIf(
597 : !stats->setList(_list),
598 : stats->errorMessage()
599 : );
600 : // Recover statistics
601 0 : Array<Double> npts, sum, sumsquared, min, max, mean, sigma;
602 0 : Array<Double> rms, fluxDensity, med, medAbsDevMed, quartile, q1, q3;
603 0 : Bool ok = true;
604 0 : auto doFlux = ! doBiweight;
605 0 : if (doFlux && this->_getImage()->imageInfo().hasMultipleBeams()) {
606 0 : if (csys.hasSpectralAxis() || csys.hasPolarizationCoordinate()) {
607 0 : Int spAxis = csys.spectralAxisNumber();
608 0 : Int poAxis = csys.polarizationAxisNumber();
609 0 : for (Int i=0; i<(Int)_axes.size(); ++i) {
610 0 : if (_axes[i] == spAxis || _axes[i] == poAxis) {
611 0 : *this->_getLog() << LogIO::WARN << "At least one cursor axis contains multiple beams. "
612 : << "You should thus use care in interpreting these statistics. Flux densities "
613 0 : << "will not be computed." << LogIO::POST;
614 0 : doFlux = false;
615 0 : break;
616 : }
617 : }
618 : }
619 : }
620 0 : if (_robust) {
621 0 : ok = stats->getStatistic(med, LatticeStatsBase::MEDIAN)
622 0 : && stats->getStatistic(
623 : medAbsDevMed, LatticeStatsBase::MEDABSDEVMED
624 : )
625 0 : && stats->getStatistic(
626 : quartile, LatticeStatsBase::QUARTILE
627 : )
628 0 : && stats->getStatistic(
629 : q1, LatticeStatsBase::Q1
630 : )
631 0 : && stats->getStatistic(
632 : q3, LatticeStatsBase::Q3
633 : );
634 : }
635 0 : ok = ok && stats->getStatistic(npts, LatticeStatsBase::NPTS)
636 0 : && stats->getStatistic(min, LatticeStatsBase::MIN)
637 0 : && stats->getStatistic(max, LatticeStatsBase::MAX)
638 0 : && stats->getStatistic(mean, LatticeStatsBase::MEAN)
639 0 : && stats->getStatistic(sigma, LatticeStatsBase::SIGMA);
640 0 : if (! doBiweight) {
641 0 : ok = ok && stats->getStatistic(sum, LatticeStatsBase::SUM)
642 0 : && stats->getStatistic(sumsquared, LatticeStatsBase::SUMSQ)
643 0 : && stats->getStatistic(rms, LatticeStatsBase::RMS);
644 : }
645 0 : ThrowIf(! ok, stats->errorMessage());
646 0 : Record statsout;
647 0 : statsout.define("npts", npts);
648 0 : statsout.define("min", min);
649 0 : statsout.define("max", max);
650 0 : statsout.define("mean", mean);
651 0 : statsout.define(SIGMA, sigma);
652 0 : if (! doBiweight) {
653 0 : statsout.define("sum", sum);
654 0 : statsout.define("sumsq", sumsquared);
655 0 : statsout.define("rms", rms);
656 : }
657 0 : if (_robust) {
658 0 : statsout.define("median", med);
659 0 : statsout.define("medabsdevmed", medAbsDevMed);
660 0 : statsout.define("quartile", quartile);
661 0 : statsout.define("q1", q1);
662 0 : statsout.define("q3", q3);
663 : }
664 0 : if (
665 : doFlux
666 0 : && stats->getStatistic(
667 : fluxDensity, LatticeStatsBase::FLUX
668 : )
669 : ) {
670 0 : statsout.define("flux", fluxDensity);
671 : }
672 0 : statsout.define("blc", blc.asVector());
673 0 : statsout.define("blcf", blcf);
674 0 : statsout.define("trc", trc.asVector());
675 0 : statsout.define("trcf", trcf);
676 0 : String tmp;
677 0 : IPosition minPos, maxPos;
678 0 : if (! doBiweight && stats->getMinMaxPos(minPos, maxPos)) {
679 0 : if (minPos.nelements() > 0) {
680 0 : statsout.define("minpos", (blc + minPos).asVector());
681 0 : tmp = CoordinateUtil::formatCoordinate(blc + minPos, csys, precis);
682 0 : statsout.define("minposf", tmp);
683 : }
684 0 : if (maxPos.nelements() > 0) {
685 0 : statsout.define("maxpos", (blc + maxPos).asVector());
686 0 : tmp = CoordinateUtil::formatCoordinate(blc + maxPos, csys, precis);
687 0 : statsout.define("maxposf", tmp);
688 : }
689 : }
690 0 : if (_list) {
691 0 : stats->showRobust(_robust);
692 0 : ThrowIf(
693 : ! stats->display(),
694 : stats->errorMessage()
695 : );
696 : }
697 0 : if (messageStore) {
698 0 : std::vector<String> messages = stats->getMessages();
699 0 : for (
700 0 : std::vector<String>::const_iterator iter=messages.begin();
701 0 : iter!=messages.end(); ++iter
702 : ) {
703 0 : messageStore->push_back(*iter + "\n");
704 : }
705 0 : stats->clearMessages();
706 : }
707 0 : return statsout;
708 : }
709 :
710 0 : template <class T> void ImageStatsCalculator<T>::_logStartup(
711 : std::vector<String> *const &messageStore, const String& myAlg,
712 : const casacore::IPosition& blc, const casacore::IPosition& trc,
713 : const casacore::String& blcf, const casacore::String trcf
714 : ) const {
715 0 : if (! _list) {
716 0 : return;
717 : }
718 0 : LogOrigin myOrigin(_class, __func__);
719 0 : *this->_getLog() << myOrigin << LogIO::NORMAL;
720 0 : String algInfo = "Statistics calculated using "
721 : + myAlg + " algorithm";
722 0 : *this->_getLog() << algInfo << LogIO::POST;
723 0 : if (messageStore) {
724 0 : messageStore->push_back(algInfo + "\n");
725 : }
726 : // Only write to the logger if the user wants it displayed.
727 0 : Vector<String> x(5);
728 0 : ostringstream y;
729 0 : x[0] = "Regions --- ";
730 0 : y << " -- bottom-left corner (pixel) [blc]: " << blc;
731 0 : x[1] = y.str();
732 0 : y.str("");
733 0 : y << " -- top-right corner (pixel) [trc]: " << trc;
734 0 : x[2] = y.str();
735 0 : y.str("");
736 0 : y << " -- bottom-left corner (world) [blcf]: " << blcf;
737 0 : x[3] = y.str();
738 0 : y.str("");
739 0 : y << " -- top-right corner (world) [trcf]: " << trcf;
740 0 : x[4] = y.str();
741 0 : for (uInt i=0; i<x.size(); ++i) {
742 0 : *this->_getLog() << x[i] << LogIO::POST;
743 0 : if (messageStore != 0) {
744 0 : messageStore->push_back(x[i] + "\n");
745 : }
746 : }
747 : }
748 :
749 755 : template <class T> void ImageStatsCalculator<T>::setRobust(Bool b) {
750 755 : _robust = b;
751 755 : }
752 :
753 0 : template <class T> casacore::Bool ImageStatsCalculator<T>::_haveRegionsChanged(
754 : ImageRegion* newRegion,
755 : ImageRegion* newMask, ImageRegion* oldRegion,
756 : ImageRegion* oldMask
757 : ) {
758 0 : Bool regionChanged = (
759 0 : newRegion != 0 && oldRegion != 0
760 0 : && (*newRegion) != (*oldRegion)
761 : )
762 0 : || (newRegion == 0 && oldRegion != 0)
763 0 : || (newRegion != 0 && oldRegion == 0
764 : );
765 0 : Bool maskChanged = (
766 0 : newMask != 0 && oldMask != 0
767 0 : && (*newMask) != (*oldMask)
768 : )
769 0 : || (newMask == 0 && oldMask != 0)
770 0 : || (newMask != 0 && oldMask == 0
771 : );
772 0 : return (regionChanged || maskChanged);
773 : }
774 :
775 : }
776 :
777 : #endif
|
