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             : #include <imageanalysis/ImageAnalysis/ImageCollapser.h>
      27             : 
      28             : #include <casacore/casa/Arrays/ArrayLogical.h>
      29             : #include <casacore/casa/BasicSL/STLIO.h>
      30             : #include <casacore/scimath/StatsFramework/ClassicalStatistics.h>
      31             : #include <casacore/scimath/Mathematics/NumericTraits.h>
      32             : #include <casacore/images/Images/ImageStatistics.h>
      33             : #include <casacore/images/Images/ImageUtilities.h>
      34             : #include <casacore/images/Images/PagedImage.h>
      35             : #include <imageanalysis/ImageAnalysis/SubImageFactory.h>
      36             : #include <casacore/images/Images/TempImage.h>
      37             : #include <casacore/lattices/Lattices/LatticeUtilities.h>
      38             : #include <casacore/lattices/LatticeMath/LatticeMathUtil.h>
      39             : 
      40             : #include <imageanalysis/ImageAnalysis/ImageMaskedPixelReplacer.h>
      41             : 
      42             : #include <memory>
      43             : 
      44             : namespace casa {
      45             : 
      46           0 : template<class T> ImageCollapser<T>::ImageCollapser(
      47             :     const casacore::String & aggString, const SPCIIT image,
      48             :     const casacore::Record * const regionRec,
      49             :     const casacore::String & maskInp, const casacore::IPosition & axes,
      50             :     casacore::Bool invertAxesSelection,
      51             :     const casacore::String & outname, casacore::Bool overwrite
      52             : ) : ImageTask<T>(
      53             :         image, "", regionRec, "", "", "",
      54             :         maskInp, outname, overwrite
      55             :     ), _invertAxesSelection(invertAxesSelection),
      56           0 :     _axes(axes), _aggType(ImageCollapserData::UNKNOWN) {
      57           0 :     _aggType = ImageCollapserData::aggregateType(aggString);
      58           0 :     this->_construct();
      59           0 :     _finishConstruction();
      60           0 : }
      61             : 
      62           0 : template<class T> ImageCollapser<T>::ImageCollapser(
      63             :     const SPCIIT image,
      64             :     const casacore::IPosition & axes, const casacore::Bool invertAxesSelection,
      65             :     const ImageCollapserData::AggregateType aggregateType,
      66             :     const casacore::String & outname, const casacore::Bool overwrite
      67             : ) : ImageTask<T>(image, "", 0, "", "", "", "", outname, overwrite),
      68             :     _invertAxesSelection(invertAxesSelection),
      69           0 :     _axes(axes), _aggType(aggregateType) {
      70           0 :     ThrowIf (
      71             :         _aggType == ImageCollapserData::UNKNOWN,
      72             :         "UNKNOWN aggregateType not allowed"
      73             :     );
      74           0 :     ThrowIf(
      75             :         ! image,
      76             :         "Cannot use a null image pointer with this constructor"
      77             :     );
      78           0 :     this->_construct();
      79           0 :     _finishConstruction();
      80           0 : }
      81             : 
      82           0 : template<class T> SPIIT ImageCollapser<T>::collapse() const {
      83           0 :     auto subImage = SubImageFactory<T>::createSubImageRO(
      84           0 :         *this->_getImage(), *this->_getRegion(), this->_getMask(),
      85             :         this->_getLog().get(), AxesSpecifier(), this->_getStretch()
      86             :     );
      87           0 :     *this->_getLog() << LogOrigin(getClass(), __func__);
      88           0 :     ThrowIf(
      89             :         ImageMask::isAllMaskFalse(*subImage),
      90             :         "All selected pixels are masked"
      91             :     );
      92           0 :     auto outCoords = subImage->coordinates();
      93           0 :     auto hasDir = outCoords.hasDirectionCoordinate();
      94           0 :     auto inShape = subImage->shape();
      95           0 :     if (_aggType == ImageCollapserData::FLUX) {
      96           0 :         _checkFlux(subImage);
      97             :     }
      98             :     // Set the compressed axis reference pixel and reference value
      99           0 :     Vector<Double> blc, trc;
     100           0 :     IPosition pixblc(inShape.nelements(), 0);
     101           0 :     auto pixtrc = inShape - 1;
     102           0 :     ThrowIf(
     103             :         ! outCoords.toWorld(blc, pixblc)
     104             :         || ! outCoords.toWorld(trc, pixtrc),
     105             :         "Could not set new coordinate values"
     106             :     );
     107           0 :     auto refValues = outCoords.referenceValue();
     108           0 :     auto refPixels = outCoords.referencePixel();
     109           0 :     auto outShape = inShape;
     110           0 :     auto degenTypes = ImageCollapserData::aggTypesSupportedDegenAxes();
     111           0 :     auto useDegenCase = degenTypes->find(_aggType) != degenTypes->end();
     112           0 :     for (const auto& axis: _axes) {
     113           0 :         refValues[axis] = (blc[axis] + trc[axis])/2;
     114           0 :         refPixels[axis] = 0;
     115           0 :         outShape[axis] = 1;
     116           0 :         useDegenCase = useDegenCase && inShape[axis] == 1;
     117             :     }
     118           0 :     ThrowIf(
     119             :         ! outCoords.setReferenceValue(refValues),
     120             :         "Unable to set reference value"
     121             :     );
     122           0 :     ThrowIf(
     123             :         ! outCoords.setReferencePixel(refPixels),
     124             :         "Unable to set reference pixel"
     125             :     );
     126           0 :     TempImage<T> tmpIm(outShape, outCoords);
     127           0 :     if (_aggType == ImageCollapserData::ZERO) {
     128           0 :         tmpIm.set(0.0);
     129             :     }
     130           0 :     else if (useDegenCase) {
     131           0 :         _doDegenerateAxesCase(tmpIm, subImage);
     132             :     }
     133           0 :     else if (
     134           0 :         _aggType == ImageCollapserData::MEDIAN
     135           0 :         || _aggType == ImageCollapserData::MADM
     136           0 :         || _aggType == ImageCollapserData::XMADM       
     137             :     ) {
     138           0 :         _doHighPerf(subImage, tmpIm);
     139             :     }
     140             :     else {
     141           0 :         _doOtherStats(tmpIm, subImage);
     142             :     }
     143           0 :     auto copied = subImage->imageInfo().hasMultipleBeams()
     144           0 :         ? _doMultipleBeams(tmpIm, subImage, hasDir, outCoords)
     145             :         : false;
     146           0 :     if (! copied) {
     147           0 :         ImageUtilities::copyMiscellaneous(tmpIm, *subImage, true);
     148             :     }
     149           0 :     if (_aggType == ImageCollapserData::FLUX) {
     150           0 :         _doFluxUnits(tmpIm, subImage);
     151             :     }
     152           0 :     return this->_prepareOutputImage(tmpIm);
     153             : }
     154             : 
     155           0 : template<class T> void ImageCollapser<T>::_checkFlux(
     156             :     SPCIIT subImage
     157             : ) const {
     158           0 :     String cant = " Cannot do flux density calculation";
     159           0 :     const auto& outCoords = subImage->coordinates();
     160           0 :     ThrowIf(
     161             :         ! outCoords.hasDirectionCoordinate(),
     162             :         "Image has no direction coordinate." + cant
     163             :     );
     164           0 :     ThrowIf(
     165             :         subImage->units().getName().contains("beam")
     166             :         && ! subImage->imageInfo().hasBeam(),
     167             :         "Image has no beam." + cant
     168             :     );
     169           0 :     auto dirAxes = outCoords.directionAxesNumbers();
     170           0 :     const auto naxes = _axes.size();
     171           0 :     for (uInt i = 0; i < naxes; ++i) {
     172           0 :         Int axis = _axes[i];
     173           0 :         ThrowIf(
     174             :             ! anyTrue(dirAxes == axis)
     175             :             && subImage->shape()[axis] > 1,
     176             :             "Specified axis " + String::toString(axis)
     177             :             + " is not a direction axis but has length > 1." + cant
     178             :         );
     179             :     }
     180           0 : }
     181             : 
     182           0 : template<class T> void ImageCollapser<T>::_doDegenerateAxesCase(
     183             :     TempImage<T>& tmpIm, SPCIIT subImage
     184             : ) const {
     185           0 :     *this->_getLog() << LogOrigin(getClass(), __func__);
     186           0 :     *this->_getLog() << LogIO::NORMAL << "All subimage axes to be "
     187             :         << "collapsed are degenerate, using algorithm optimized for "
     188           0 :         << "that case." << LogIO::POST;
     189           0 :     ThrowIf(
     190             :         _aggType == ImageCollapserData::STDDEV
     191             :         || _aggType == ImageCollapserData::VARIANCE,
     192             :         "Cannot compute "
     193             :         + ImageCollapserData::funcNameMap()->find(_aggType)->second
     194             :         + " for single sample data sets"
     195             :     );
     196           0 :     if (
     197           0 :         _aggType == ImageCollapserData::MAX
     198           0 :         || _aggType == ImageCollapserData::MEAN
     199           0 :         || _aggType == ImageCollapserData::MEDIAN
     200           0 :         || _aggType == ImageCollapserData::MIN
     201           0 :         || _aggType == ImageCollapserData::SUM
     202             :     ) {
     203             :         // Straight copy
     204           0 :         this->_copyData(tmpIm, *subImage);
     205             :     }
     206           0 :     else if (_aggType == ImageCollapserData::NPTS) {
     207           0 :         tmpIm.set(1.0);
     208             :     }
     209           0 :     else if (
     210           0 :         _aggType == ImageCollapserData::MADM
     211           0 :         || _aggType == ImageCollapserData::XMADM
     212             :     ) {
     213           0 :         tmpIm.set(0.0);
     214             :     }
     215           0 :     else if (_aggType == ImageCollapserData::RMS) {
     216           0 :         this->_copyData(tmpIm, LatticeExpr<T>(abs(*subImage)));
     217             :     }
     218             :     else {
     219           0 :         ThrowCc(
     220             :             "Logic error: "
     221             :             + ImageCollapserData::funcNameMap()->find(_aggType)->second
     222             :             + " erroneously not supported for degenerate axis case. Please "
     223             :             + "file a bug report and include this message"
     224             :         );
     225             :     }
     226           0 :     if (subImage->isMasked() && ! ImageMask::isAllMaskTrue(*subImage)) {
     227           0 :         if (! tmpIm.isMasked()) {
     228           0 :             TempLattice<Bool> mask(tmpIm.shape());
     229           0 :             this->_copyMask(mask, *subImage);
     230           0 :             tmpIm.attachMask(mask);
     231             :         }
     232             :         // This works because the underlying pixel data are cloned by reference
     233           0 :         SPIIT myclone(tmpIm.cloneII());
     234           0 :         ImageMaskedPixelReplacer<T> impr(myclone);
     235           0 :         impr.replace("0", False, False);
     236             :     }
     237           0 : }
     238             : 
     239           0 : template<class T> void ImageCollapser<T>::_doFluxUnits(
     240             :     TempImage<T>& tmpIm, const std::shared_ptr<const SubImage<T>> subImage
     241             : ) const {
     242             :     // get the flux units right
     243           0 :      auto sbunit = subImage->units().getName();
     244           0 :      String unit;
     245           0 :      if (sbunit.contains("K")) {
     246           0 :          casacore::String areaUnit = "arcsec2";
     247           0 :          unit = sbunit + "." + areaUnit;
     248             :      }
     249             :      else {
     250           0 :          unit = "Jy";
     251           0 :          if (sbunit.contains("/beam")) {
     252           0 :              uInt iBeam = sbunit.find("/beam");
     253           0 :              unit = sbunit.substr(0, iBeam) + sbunit.substr(iBeam + 5);
     254             :          }
     255             :      }
     256           0 :      tmpIm.setUnits(unit);
     257           0 : }
     258             : 
     259           0 : template<class T> void ImageCollapser<T>::_doHighPerf(
     260             :     SPCIIT image, casacore::TempImage<T>& outImage
     261             : ) const {
     262           0 :     auto doMedian = _aggType == ImageCollapserData::MEDIAN;
     263           0 :     auto doMADM = _aggType == ImageCollapserData::MADM
     264           0 :         || _aggType == ImageCollapserData::XMADM;
     265           0 :     ThrowIf(
     266             :         ! doMedian && ! doMADM,
     267             :         "Logic error, unsupported aggregate type "
     268             :         + String(ImageCollapserData::funcNameMap()->at((uInt)_aggType)) + " for method "
     269             :         + String(__func__)
     270             :     );
     271           0 :     IPosition cursorShape(image->ndim(), 1);
     272           0 :     for (uInt i = 0; i < cursorShape.size(); ++i) {
     273           0 :         for (uInt j = 0; j < _axes.size(); ++j) {
     274           0 :             if (_axes[j] == i) {
     275           0 :                 cursorShape[i] = image->shape()[i];
     276           0 :                 break;
     277             :             }
     278             :         }
     279             :     }
     280           0 :     LatticeStepper stepper(image->shape(), cursorShape);
     281           0 :     std::unique_ptr<Array<Bool>> outMask;
     282             :     // accumtype being the same precision as the input data type is ok here,
     283             :     // since we are only computing the median/madm and not actually accumulating
     284             :     ClassicalStatistics<
     285             :         T, typename Array<T>::const_iterator, Array<Bool>::const_iterator
     286           0 :     > stats;
     287           0 :     auto hasMaskedPixels = ! ImageMask::isAllMaskTrue(*image);
     288           0 :     for (stepper.reset(); !stepper.atEnd(); stepper++) {
     289           0 :         Slicer slicer(
     290             :             stepper.position(), stepper.endPosition(), casacore::Slicer::endIsLast
     291             :         );
     292           0 :         auto data = image->getSlice(slicer);
     293           0 :         Bool isMasked = False;
     294           0 :         Array<Bool> maskSlice;
     295           0 :         if (hasMaskedPixels) {
     296           0 :             maskSlice = image->getMaskSlice(slicer);
     297           0 :             isMasked = ! allTrue(maskSlice);
     298             :         }
     299           0 :         if (isMasked) {
     300           0 :             if (! anyTrue(maskSlice)) {
     301           0 :                 if (! outMask) {
     302           0 :                     outMask.reset(new Array<Bool>(outImage.shape(), true));
     303             :                 }
     304           0 :                 (*outMask)(stepper.position()) = false;
     305           0 :                 outImage.putAt(0, stepper.position());
     306             :             }
     307           0 :             else if (! allTrue(maskSlice)) {
     308           0 :                 stats.setData(data.begin(), maskSlice.begin(), data.size());
     309           0 :                 if (doMedian) {
     310           0 :                     outImage.putAt(stats.getMedian(), stepper.position());
     311             :                 }
     312           0 :                 else if (doMADM) {
     313           0 :                     auto x = stats.getMedianAbsDevMed();
     314           0 :                     if (_aggType == ImageCollapserData::XMADM) {
     315           0 :                         x *= C::probit_3_4;
     316             :                     }
     317           0 :                     outImage.putAt(x, stepper.position());
     318             :                 }
     319             :             }
     320             :         }
     321             :         else {
     322           0 :             stats.setData(data.begin(), data.size());
     323           0 :             if (doMedian) {
     324           0 :                 outImage.putAt(stats.getMedian(), stepper.position());
     325             :             }
     326           0 :             else if (doMADM) {
     327           0 :                 auto x = stats.getMedianAbsDevMed();
     328           0 :                 if (_aggType == ImageCollapserData::XMADM) {
     329           0 :                     x *= C::probit_3_4;
     330             :                 }
     331           0 :                 outImage.putAt(x, stepper.position());
     332             :             }
     333             :         }
     334             :     }
     335           0 :     if (outMask) {
     336           0 :         outImage.attachMask(ArrayLattice<Bool>(*outMask));
     337             :     }
     338           0 : }
     339             : 
     340           0 : template<class T> Bool ImageCollapser<T>::_doMultipleBeams(
     341             :     TempImage<T>& tmpIm, SPCIIT subImage, Bool hasDir,
     342             :     const CoordinateSystem & outCoords
     343             : ) const {
     344           0 :     auto naxes = _axes.size();
     345           0 :     auto dirAxesOnlyCollapse = hasDir && naxes == 2;
     346           0 :     if (dirAxesOnlyCollapse) {
     347           0 :         auto dirAxes = outCoords.directionAxesNumbers();
     348           0 :         dirAxesOnlyCollapse = (_axes[0] == dirAxes[0] && _axes[1] == dirAxes[1])
     349           0 :                               || (_axes[1] == dirAxes[0] && _axes[0] == dirAxes[1]);
     350             :     }
     351           0 :     if (! dirAxesOnlyCollapse) {
     352             :         // check for degeneracy of spectral or polarization axes
     353           0 :         auto specAxis = outCoords.spectralAxisNumber(false);
     354           0 :         auto polAxis = outCoords.polarizationAxisNumber(false);
     355           0 :         dirAxesOnlyCollapse = true;
     356           0 :         auto shape = subImage->shape();
     357           0 :         for (uInt i = 0; i < naxes; ++i) {
     358           0 :             auto axis = _axes[i];
     359           0 :             if (
     360           0 :                 (axis == specAxis || axis == polAxis)
     361           0 :                 && shape[axis] > 1
     362             :             ) {
     363           0 :                 dirAxesOnlyCollapse = false;
     364           0 :                 break;
     365             :             }
     366             :         }
     367             :     }
     368           0 :     if (! dirAxesOnlyCollapse) {
     369           0 :         LogOrigin lor(getClass(), __func__);
     370           0 :         String msg = "Input image has per plane beams "
     371             :             "but the collapse is not done exclusively along the direction axes. "
     372             :             "The output image will arbitrarily have a single beam which "
     373             :             "is the first beam available in the subimage."
     374             :             "Thus, the image planes will not be convolved to a common "
     375             :             "restoring beam before collapsing. If, however, this is desired, "
     376             :             "then run the task imsmooth or the tool method ia.convolve2d() first, "
     377             :             "and use the output image of that as the input for collapsing.";
     378           0 :         *this->_getLog() << lor << LogIO::WARN << msg << LogIO::POST;
     379           0 :         this->addHistory(lor, msg);
     380           0 :         ImageUtilities::copyMiscellaneous(tmpIm, *subImage, false);
     381           0 :         auto info = subImage->imageInfo();
     382           0 :         auto beam = *(info.getBeamSet().getBeams().begin());
     383           0 :         info.removeRestoringBeam();
     384           0 :         info.setRestoringBeam(beam);
     385           0 :         tmpIm.setImageInfo(info);
     386           0 :         return true;
     387             :     }
     388           0 :     return false;
     389             : }
     390             : 
     391           0 : template<class T> void ImageCollapser<T>::_doOtherStats(
     392             :     TempImage<T>& tmpIm, SPCIIT subImage
     393             : ) const {
     394           0 :     T npixPerBeam = 1;
     395           0 :     if (_aggType == ImageCollapserData::SQRTSUM_NPIX_BEAM) {
     396           0 :         const auto& info = subImage->imageInfo();
     397           0 :         if (! info.hasBeam()) {
     398           0 :             *this->_getLog() << casacore::LogIO::WARN
     399             :                 << "Image has no beam, will use sqrtsum method"
     400           0 :                 << casacore::LogIO::POST;
     401             :         }
     402           0 :         else if (info.hasMultipleBeams()) {
     403           0 :             *this->_getLog() << casacore::LogIO::WARN
     404             :                 << "Function sqrtsum_npix_beam does not support multiple beams, will"
     405             :                 << "use sqrtsum method instead"
     406           0 :                 << casacore::LogIO::POST;
     407             :         }
     408             :         else {
     409           0 :             npixPerBeam = info.getBeamAreaInPixels(
     410           0 :                 -1, -1, subImage->coordinates().directionCoordinate()
     411             :             );
     412             :         }
     413             :     }
     414           0 :     _doLowPerf(tmpIm, subImage, npixPerBeam);
     415           0 : }
     416             : 
     417           0 : template<class T> void ImageCollapser<T>::_doLowPerf(
     418             :     TempImage<T>& tmpIm, SPCIIT subImage, T npixPerBeam
     419             : ) const {
     420             :     // flux or mask with one or more false values, must use lower performance methods
     421           0 :     auto lattStatType = _getStatsType();
     422           0 :     Array<T> data;
     423           0 :     Array<Bool> mask;
     424           0 :     if (_aggType == ImageCollapserData::FLUX) {
     425           0 :         ImageStatistics<T> stats(*subImage, false);
     426           0 :         stats.setAxes(_axes.asVector());
     427           0 :         if (
     428           0 :             ! stats.getConvertedStatistic(data, lattStatType, false)
     429             :         ) {
     430           0 :             ostringstream oss;
     431           0 :             oss << "Unable to calculate flux density: "
     432           0 :                 << stats.getMessages();
     433           0 :             ThrowCc(oss.str());
     434             :         }
     435           0 :         mask.resize(data.shape());
     436           0 :         mask.set(true);
     437             :     }
     438             :     else {
     439           0 :         LatticeMathUtil::collapse(
     440           0 :             data, mask, _axes, *subImage, false,
     441             :             true, true, lattStatType
     442             :         );
     443           0 :         if (
     444           0 :             _aggType == ImageCollapserData::SQRTSUM
     445           0 :             || _aggType == ImageCollapserData::SQRTSUM_NPIX
     446           0 :             || _aggType == ImageCollapserData::SQRTSUM_NPIX_BEAM
     447             :         ) {
     448           0 :             _zeroNegatives(data);
     449           0 :             data = sqrt(data);
     450           0 :             if (_aggType == ImageCollapserData::SQRTSUM_NPIX) {
     451           0 :                 auto npts = data.copy();
     452           0 :                 LatticeMathUtil::collapse(
     453           0 :                     npts, mask, _axes, *subImage, false,
     454             :                     true, true, LatticeStatsBase::NPTS
     455             :                 );
     456           0 :                 data /= npts;
     457             :             }
     458           0 :             else if (_aggType == ImageCollapserData::SQRTSUM_NPIX_BEAM) {
     459           0 :                 data /= npixPerBeam;
     460             :             }
     461             :         }
     462             :     }
     463           0 :     auto dataCopy = (_axes.size() <= 1)
     464           0 :         ? data : data.addDegenerate(_axes.size() - 1);
     465           0 :     IPosition newOrder(tmpIm.ndim(), -1);
     466           0 :     auto nAltered = _axes.size();
     467           0 :     auto nUnaltered = tmpIm.ndim() - nAltered;
     468           0 :     auto alteredCount = nUnaltered;
     469           0 :     auto unAlteredCount = 0;
     470           0 :     const auto ndim = tmpIm.ndim();
     471           0 :     const auto naxes = nAltered;
     472           0 :     for (uInt i = 0; i < ndim; ++i) {
     473           0 :         for (uInt j = 0; j < naxes; ++j) {
     474           0 :             if (i == _axes[j]) {
     475           0 :                 newOrder[i] = alteredCount;
     476           0 :                 alteredCount++;
     477           0 :                 break;
     478             :             }
     479             :         }
     480           0 :         if (newOrder[i] < 0) {
     481           0 :             newOrder[i] = unAlteredCount;
     482           0 :             ++unAlteredCount;
     483             :         }
     484             :     }
     485           0 :     tmpIm.put(reorderArray(dataCopy, newOrder));
     486           0 :     if (! allTrue(mask)) {
     487           0 :         auto maskCopy = (
     488           0 :             _axes.size() <= 1) ? mask
     489           0 :                 : mask.addDegenerate(_axes.size() - 1
     490             :         );
     491           0 :         auto mCopy = reorderArray(maskCopy, newOrder);
     492           0 :         tmpIm.attachMask(ArrayLattice<Bool>(mCopy));
     493             :     }
     494           0 : }
     495             : 
     496           0 : template<class T> LatticeStatsBase::StatisticsTypes ImageCollapser<T>::_getStatsType() const {
     497           0 :     auto lattStatType = LatticeStatsBase::NACCUM;
     498           0 :     switch (_aggType) {
     499           0 :     case ImageCollapserData::FLUX:
     500           0 :         lattStatType = LatticeStatsBase::FLUX;
     501           0 :         break;
     502           0 :     case ImageCollapserData::MAX:
     503           0 :         lattStatType = LatticeStatsBase::MAX;
     504           0 :         break;
     505           0 :     case ImageCollapserData::MEAN:
     506           0 :         lattStatType = LatticeStatsBase::MEAN;
     507           0 :         break;
     508           0 :     case ImageCollapserData::MIN:
     509           0 :         lattStatType = LatticeStatsBase::MIN;
     510           0 :         break;
     511           0 :     case ImageCollapserData::NPTS:
     512           0 :         lattStatType = LatticeStatsBase::NPTS;
     513           0 :         break;
     514           0 :     case ImageCollapserData::RMS:
     515           0 :         lattStatType = LatticeStatsBase::RMS;
     516           0 :         break;
     517           0 :     case ImageCollapserData::STDDEV:
     518           0 :         lattStatType = LatticeStatsBase::SIGMA;
     519           0 :         break;
     520           0 :     case ImageCollapserData::SQRTSUM:
     521             :     case ImageCollapserData::SQRTSUM_NPIX:
     522             :     case ImageCollapserData::SQRTSUM_NPIX_BEAM:
     523             :     case ImageCollapserData::SUM:
     524           0 :         lattStatType = LatticeStatsBase::SUM;
     525           0 :         break;
     526           0 :     case ImageCollapserData::VARIANCE:
     527           0 :         lattStatType = LatticeStatsBase::VARIANCE;
     528           0 :         break;
     529           0 :     case ImageCollapserData::MEDIAN:
     530             :     case ImageCollapserData::ZERO:
     531             :     case ImageCollapserData::UNKNOWN:
     532             :     default:
     533           0 :         ThrowCc(
     534             :             "Logic error. Should never have gotten the the bottom of the switch statement"
     535             :         );
     536             :         break;
     537             :     }
     538           0 :     return lattStatType;
     539             : }
     540             : 
     541           0 : template<class T> void ImageCollapser<T>::_zeroNegatives(Array<T>& arr) {
     542           0 :     auto iter = arr.begin();
     543           0 :     if (isComplex(whatType<T>()) || allGE(arr, (T)0)) {
     544           0 :         return;
     545             :     }
     546           0 :     auto end = arr.end();
     547           0 :     for (; iter != end; ++iter) {
     548           0 :         if (*iter < 0) {
     549           0 :             *iter = 0;
     550             :         }
     551             :     }
     552             : }
     553             : 
     554           0 : template<class T> void ImageCollapser<T>::_finishConstruction() {
     555           0 :     for (
     556           0 :         casacore::IPosition::const_iterator iter = _axes.begin();
     557           0 :         iter != _axes.end(); iter++
     558             :     ) {
     559           0 :         ThrowIf(
     560             :             *iter >= this->_getImage()->ndim(),
     561             :             "Specified zero-based axis (" + casacore::String::toString(*iter)
     562             :             + ") must be less than the number of axes in " + this->_getImage()->name()
     563             :             + "(" + casacore::String::toString(this->_getImage()->ndim()) + ")"
     564             :         );
     565             :     }
     566           0 :     _invert();
     567           0 : }
     568             : 
     569           0 : template<class T> void ImageCollapser<T>::_invert() {
     570           0 :     if (_invertAxesSelection) {
     571           0 :         casacore::IPosition x = casacore::IPosition::otherAxes(this->_getImage()->ndim(), _axes);
     572           0 :         _axes.resize(x.size());
     573           0 :         _axes = x;
     574             :     }
     575           0 : }
     576             : 
     577             : }