Line data    Source code

       1             : #include <imageanalysis/ImageAnalysis/StatImageCreator.h>
       2             : 
       3             : #include <imageanalysis/Annotations/AnnCenterBox.h>
       4             : #include <imageanalysis/Annotations/AnnCircle.h>
       5             : #include <casacore/lattices/LEL/LatticeExpr.h>
       6             : 
       7             : namespace casa {
       8             : 
       9          34 : StatImageCreator::StatImageCreator(
      10             :     const SPCIIF image, const Record *const region,
      11             :         const String& mask, const String& outname, Bool overwrite
      12          34 : ) : ImageStatsBase(image, region, mask, outname, overwrite) {
      13          34 :     this->_construct();
      14          34 :         auto da = _getImage()->coordinates().directionAxesNumbers();
      15          34 :     _dirAxes[0] = da[0];
      16          34 :     _dirAxes[1] = da[1];
      17          34 :         setAnchorPosition(0, 0);
      18          34 : }
      19             : 
      20           1 : void StatImageCreator::setRadius(const Quantity& radius) {
      21           1 :     ThrowIf(
      22             :         ! (radius.getUnit().startsWith("pix") || radius.isConform("rad")),
      23             :         "radius units " + radius.getUnit() + " must be pixel or angular"
      24             :     );
      25           1 :     _xlen = radius;
      26           1 :     _ylen = Quantity(0,"");
      27           1 : }
      28             : 
      29          33 : void StatImageCreator::setRectangle(
      30             :     const Quantity& xLength, const Quantity& yLength
      31             : ) {
      32          33 :     ThrowIf(
      33             :         ! (xLength.getUnit().startsWith("pix") || xLength.isConform("rad")),
      34             :         "xLength units " + xLength.getUnit() + " must be pixel or angular"
      35             :     );
      36          33 :     ThrowIf(
      37             :         ! (yLength.getUnit().startsWith("pix") || yLength.isConform("rad")),
      38             :         "xLength units " + yLength.getUnit() + " must be pixel or angular"
      39             :     );
      40          33 :     _xlen = xLength;
      41          33 :     _ylen = yLength;
      42          33 : }
      43             : 
      44          68 : void StatImageCreator::setAnchorPosition(Int x, Int y) {
      45          68 :     Vector<Double> anchorPixel(_getImage()->ndim(), 0);
      46          68 :     anchorPixel[_dirAxes[0]] = x;
      47          68 :     anchorPixel[_dirAxes[1]] = y;
      48          68 :     _anchor = _getImage()->coordinates().toWorld(anchorPixel);
      49          68 : }
      50             : 
      51          30 : void StatImageCreator::useReferencePixelAsAnchor() {
      52          30 :     const auto refPix = _getImage()->coordinates().referencePixel();
      53          30 :     Int x = round(refPix[_dirAxes[0]]);
      54          30 :     Int y = round(refPix[_dirAxes[1]]);
      55          90 :     *_getLog() << LogIO::NORMAL << LogOrigin("StatImageCreator", __func__)
      56             :         << "Anchor being set at pixel [" << x << "," << y
      57          60 :         << "], at/near image reference pixel." << LogIO::POST;
      58          30 :     setAnchorPosition(x, y);
      59          30 : }
      60             : 
      61          34 : void StatImageCreator::setGridSpacing(uInt x, uInt y) {
      62          34 :     _grid.first = x;
      63          34 :     _grid.second = y;
      64          34 : }
      65             : 
      66          34 : SPIIF StatImageCreator::compute() {
      67          34 :     static const AxesSpecifier dummyAxesSpec;
      68          34 :     *_getLog() << LogOrigin(getClass(), __func__);
      69          34 :     auto mylog = _getLog().get();
      70             :     auto subImageRO = SubImageFactory<Float>::createSubImageRO(
      71          34 :         *_getImage(), *_getRegion(), _getMask(),
      72          34 :         mylog, dummyAxesSpec, _getStretch()
      73          68 :     );
      74             :     auto subImageRW = SubImageFactory<Float>::createImage(
      75         102 :         *_getImage(), "", *_getRegion(), _getMask(),
      76          34 :         dummyAxesSpec, False, False, _getStretch()
      77         102 :     );
      78          34 :     _doMask = ! ImageMask::isAllMaskTrue(*subImageRO);
      79          34 :     const auto imshape = subImageRO->shape();
      80          34 :     const auto xshape = imshape[_dirAxes[0]];
      81          34 :     const auto yshape = imshape[_dirAxes[1]];
      82          34 :     const auto& csys = subImageRO->coordinates();
      83          34 :     auto anchorPixel = csys.toPixel(_anchor);
      84          34 :     Int xanchor = rint(anchorPixel[_dirAxes[0]]);
      85          34 :     Int yanchor = rint(anchorPixel[_dirAxes[1]]);
      86             :     // ensure xanchor and yanchor are positive
      87          34 :     if (xanchor < 0) {
      88             :         // ugh, mod of a negative number in C++ doesn't do what I want it to
      89             :         // integer division
      90           0 :         xanchor += (abs(xanchor)/_grid.first + 1)*_grid.first;
      91             :     }
      92          34 :     if (yanchor < 0) {
      93           0 :         yanchor += (abs(yanchor)/_grid.second + 1)*_grid.second;
      94             :     }
      95             :     // xstart and ystart are the pixel location in the
      96             :     // subimage of the lower left corner of the grid,
      97             :     // ie they are the pixel location of the grid point
      98             :     // with the smallest non-negative x and y values
      99          34 :     Int xstart = xanchor % _grid.first;
     100          34 :     Int ystart = yanchor % _grid.second;
     101          34 :     if (xstart < 0) {
     102           0 :         xstart += _grid.first;
     103             :     }
     104          34 :     if (ystart < 0) {
     105           0 :         ystart += _grid.second;
     106             :     }
     107          34 :     Array<Float> tmp;
     108             :     // integer division
     109          34 :     auto nxpts = (xshape - xstart)/_grid.first;
     110          34 :     if ((xshape - xstart) % _grid.first != 0) {
     111           4 :         ++nxpts;
     112             :     }
     113          34 :     auto nypts = (yshape - ystart)/ _grid.second;
     114          34 :     if ((yshape - ystart) % _grid.second != 0) {
     115           3 :         ++nypts;
     116             :     }
     117          34 :     IPosition storeShape = imshape;
     118          34 :     storeShape[_dirAxes[0]] = nxpts;
     119          34 :     storeShape[_dirAxes[1]] = nypts;
     120          34 :     auto interpolate = _grid.first > 1 || _grid.second > 1;
     121          34 :     TempImage<Float>* writeTo = dynamic_cast<TempImage<Float> *>(subImageRW.get());
     122          34 :     std::unique_ptr<TempImage<Float>> store;
     123          34 :     if (interpolate) {
     124           6 :         store.reset(new TempImage<Float>(storeShape, csys));
     125           6 :         store->set(0.0);
     126           6 :         if (_doMask) {
     127           0 :             store->attachMask(ArrayLattice<Bool>(storeShape));
     128           0 :             store->pixelMask().set(True);
     129             :         }
     130           6 :         writeTo = store.get();
     131             :     }
     132          34 :     _computeStat(*writeTo, subImageRO, nxpts, nypts, xstart, ystart);
     133          34 :     if (_doProbit) {
     134           0 :         writeTo->copyData((LatticeExpr<Float>)(*writeTo * C::probit_3_4));
     135             :     }
     136          34 :     if (interpolate) {
     137          12 :         _doInterpolation(
     138           6 :             subImageRW, *store, subImageRO,
     139             :             nxpts, nypts, xstart,  ystart
     140             :         );
     141             :     }
     142          34 :     auto res = _prepareOutputImage(*subImageRW);
     143          68 :     return res;
     144          34 : }
     145             : 
     146          34 : void StatImageCreator::_computeStat(
     147             :     TempImage<Float>& writeTo, SPCIIF subImage,
     148             :     uInt nxpts, uInt nypts, Int xstart, Int ystart
     149             : ) {
     150          34 :     std::shared_ptr<Array<Bool>> regionMask;
     151          34 :     uInt xBlcOff = 0;
     152          34 :     uInt yBlcOff = 0;
     153          34 :     uInt xChunkSize = 0;
     154          34 :     uInt yChunkSize = 0;
     155          34 :     _nominalChunkInfo(
     156             :         regionMask,  xBlcOff, yBlcOff,
     157             :         xChunkSize, yChunkSize, subImage
     158             :     );
     159          34 :     Array<Bool> regMaskCopy;
     160          34 :     if (regionMask) {
     161           1 :         regMaskCopy = regionMask->copy();
     162           1 :         if (! writeTo.hasPixelMask()) {
     163           1 :             ArrayLattice<Bool> mymask(writeTo.shape());
     164           1 :             mymask.set(True);
     165           1 :             writeTo.attachMask(mymask);
     166           1 :         }
     167             :     }
     168          34 :     auto imshape = subImage->shape();
     169          34 :     auto ndim = imshape.size();
     170          34 :     IPosition chunkShape(ndim, 1);
     171          34 :     chunkShape[_dirAxes[0]] = xChunkSize;
     172          34 :     chunkShape[_dirAxes[1]] = yChunkSize;
     173          34 :     auto xsize = imshape[_dirAxes[0]];
     174          34 :     auto ysize = imshape[_dirAxes[1]];
     175          34 :     IPosition planeShape(imshape.size(), 1);
     176          34 :     planeShape[_dirAxes[0]] = xsize;
     177          34 :     planeShape[_dirAxes[1]] = ysize;
     178             :     RO_MaskedLatticeIterator<Float> lattIter (
     179          34 :         *subImage, planeShape, True
     180          34 :     );
     181          34 :     String algName;
     182          34 :     auto alg = _getStatsAlgorithm(algName);
     183          34 :     std::set<StatisticsData::STATS> statToCompute;
     184          34 :     statToCompute.insert(_statType);
     185          34 :     alg->setStatsToCalculate(statToCompute);
     186          34 :     auto ngrid = nxpts*nypts;
     187          34 :     auto nPts = ngrid;
     188          34 :     IPosition loopAxes;
     189         165 :     for (uInt i=0; i<ndim; ++i) {
     190         131 :         if (i != _dirAxes[0] && i != _dirAxes[1]) {
     191          63 :             loopAxes.append(IPosition(1, i));
     192          63 :             nPts *= imshape[i];
     193             :         }
     194             :     }
     195          34 :     auto doCircle = _ylen.getValue() <= 0;
     196          68 :     *_getLog() << LogOrigin(getClass(), __func__) << LogIO::NORMAL
     197          34 :         << "Using ";
     198          34 :     if (doCircle) {
     199           1 :         *_getLog() << "circular region of radius " << _xlen;
     200             :     }
     201             :     else {
     202          66 :         *_getLog() << "rectangular region of specified dimensions " << _xlen
     203          33 :             << " x " << _ylen;
     204             :     }
     205          68 :     *_getLog() << " (because of centering and "
     206             :             << "rounding to use whole pixels, actual dimensions of bounding box are "
     207          34 :             << xChunkSize << " pix x " << yChunkSize << " pix";
     208          34 :     if (doCircle) {
     209           2 :         *_getLog() << " and there are " << ntrue(*regionMask)
     210           1 :             << " good pixels in the circle that are being used";
     211             :     }
     212          68 :     *_getLog() << ") to choose pixels for computing " << _statName
     213             :         << " using the " << algName << " algorithm around each of "
     214          34 :         << ngrid << " grid points in " << (nPts/ngrid) << " planes." << LogIO::POST;
     215          34 :     auto imageChunkShape = chunkShape;
     216          34 :     _doStatsLoop(
     217             :         writeTo, lattIter, nxpts, nypts, xstart, ystart,
     218             :         xBlcOff, yBlcOff, xChunkSize, yChunkSize, imshape,
     219             :         chunkShape, regionMask, alg, regMaskCopy, loopAxes, nPts
     220             :     );
     221          34 : }
     222             : 
     223          34 : void StatImageCreator::_doStatsLoop(
     224             :     TempImage<Float>& writeTo, RO_MaskedLatticeIterator<Float>& lattIter,
     225             :     uInt nxpts, uInt nypts, Int xstart, Int ystart, uInt xBlcOff,
     226             :     uInt yBlcOff, uInt xChunkSize, uInt yChunkSize, const IPosition& imshape,
     227             :     const IPosition& chunkShape, std::shared_ptr<Array<Bool>> regionMask,
     228             :     std::shared_ptr<
     229             :         StatisticsAlgorithm<
     230             :             Double, Array<Float>::const_iterator, Array<Bool>::const_iterator,
     231             :             Array<Float>::const_iterator
     232             :         >
     233             :     >& alg, const Array<Bool>& regMaskCopy, const IPosition& loopAxes, uInt nPts
     234             : ) {
     235          34 :     auto ximshape = imshape[_dirAxes[0]];
     236          34 :     auto yimshape = imshape[_dirAxes[1]];
     237          34 :     auto ndim = imshape.size();
     238          34 :     IPosition planeBlc(ndim, 0);
     239          34 :     auto& xPlaneBlc = planeBlc[_dirAxes[0]];
     240          34 :     auto& yPlaneBlc = planeBlc[_dirAxes[1]];
     241          34 :     auto planeTrc = planeBlc + chunkShape - 1;
     242          34 :     auto& xPlaneTrc = planeTrc[_dirAxes[0]];
     243          34 :     auto& yPlaneTrc = planeTrc[_dirAxes[1]];
     244          34 :     uInt nCount = 0;
     245          34 :     auto hasLoopAxes = ! loopAxes.empty();
     246          34 :     ProgressMeter pm(0, nPts, "Processing stats at grid points");
     247             :     // Vector rather than IPosition because blc values can be negative
     248          34 :     Vector<Int> blc(ndim, 0);
     249          34 :     blc[_dirAxes[0]] = xstart - xBlcOff;
     250          34 :     blc[_dirAxes[1]] = ystart - yBlcOff;
     251         209 :     for (lattIter.atStart(); ! lattIter.atEnd(); ++lattIter) {
     252         175 :         auto plane = lattIter.cursor();
     253         175 :         auto lattMask = lattIter.getMask();
     254         175 :         auto checkGrid = ! allTrue(lattMask);
     255         175 :         auto outPos = lattIter.position();
     256         175 :         auto& xOutPos = outPos[_dirAxes[0]];
     257         175 :         auto& yOutPos = outPos[_dirAxes[1]];
     258             :         // inPos is the position of the current grid point
     259             :         // This is the position in the current chunk, not the entire lattice
     260         175 :         auto inPos = plane.shape() - 1;
     261         175 :         inPos[_dirAxes[0]] = xstart;
     262         175 :         inPos[_dirAxes[1]] = ystart;
     263         175 :         auto& xInPos = inPos[_dirAxes[0]];
     264         175 :         auto& yInPos = inPos[_dirAxes[1]];
     265         175 :         Int yblc = ystart - yBlcOff;
     266       42386 :         for (uInt yCount=0; yCount<nypts; ++yCount, yblc+=_grid.second, yInPos+=_grid.second) {
     267       42211 :             yOutPos = yCount;
     268       42211 :             yPlaneBlc = max(0, yblc);
     269       84422 :             yPlaneTrc = min(
     270       42211 :                 yblc + (Int)yChunkSize - 1, (Int)yimshape - 1
     271             :             );
     272       42211 :             IPosition regMaskStart(ndim, 0);
     273       42211 :             auto& xRegMaskStart = regMaskStart[_dirAxes[0]];
     274       42211 :             auto& yRegMaskStart = regMaskStart[_dirAxes[1]];
     275       42211 :             auto regMaskLength = regMaskStart;
     276       42211 :             auto& xRegMaskLength = regMaskLength[_dirAxes[0]];
     277       42211 :             auto& yRegMaskLength = regMaskLength[_dirAxes[1]];
     278       42211 :             Bool yDoMaskSlice = False;
     279       42211 :             if (regionMask) {
     280           5 :                 regMaskLength = regionMask->shape();
     281           5 :                 if (yblc < 0) {
     282           1 :                     yRegMaskStart = -yblc;
     283           1 :                     yRegMaskLength += yblc;
     284           1 :                     yDoMaskSlice = True;
     285             :                 }
     286           4 :                 else if (yblc + yChunkSize > yimshape) {
     287           1 :                     yRegMaskLength = yimshape - yblc;
     288           1 :                     yDoMaskSlice = True;
     289             :                 }
     290             :             }
     291       42211 :             xInPos = xstart;
     292       42211 :             Int xblc = xstart - xBlcOff;
     293    10933570 :             for (uInt xCount=0; xCount<nxpts; ++xCount, xblc+=_grid.first, xInPos+=_grid.first) {
     294    10891359 :                 Float res = 0;
     295    10891359 :                 xOutPos = xCount;
     296    10891359 :                 if (checkGrid && ! lattMask(inPos)) {
     297             :                     // grid point is masked, no computation should be done
     298     3014675 :                     writeTo.putAt(res, outPos);
     299     3014675 :                     writeTo.pixelMask().putAt(False, outPos);
     300             :                 }
     301             :                 else {
     302     7876684 :                     xPlaneBlc = max(0, xblc);
     303    15753368 :                     xPlaneTrc = min(
     304     7876684 :                         xblc + (Int)xChunkSize - 1, (Int)ximshape - 1
     305             :                     );
     306     7876684 :                     std::shared_ptr<Array<Bool>> subRegionMask;
     307     7876684 :                     if (regionMask) {
     308          25 :                         auto doMaskSlice = yDoMaskSlice;
     309          25 :                         xRegMaskStart = 0;
     310          25 :                         if (xblc < 0) {
     311           5 :                             xRegMaskStart = -xblc;
     312           5 :                             xRegMaskLength = regionMask->shape()[_dirAxes[0]] + xblc;
     313           5 :                             doMaskSlice = True;
     314             :                         }
     315          20 :                         else if (xblc + xChunkSize > ximshape) {
     316           5 :                             regMaskLength[_dirAxes[0]] = ximshape - xblc;
     317           5 :                             doMaskSlice = True;
     318             :                         }
     319             :                         else {
     320          15 :                             xRegMaskLength = xChunkSize;
     321             :                         }
     322          25 :                         if (doMaskSlice) {
     323          16 :                             Slicer sl(regMaskStart, regMaskLength);
     324          16 :                             subRegionMask.reset(new Array<Bool>(regMaskCopy(sl)));
     325          16 :                         }
     326             :                         else {
     327           9 :                             subRegionMask.reset(new Array<Bool>(regMaskCopy));
     328             :                         }
     329             :                     }
     330     7876684 :                     auto maskChunk = lattMask(planeBlc, planeTrc).copy();
     331     7876684 :                     if (subRegionMask) {
     332          25 :                         maskChunk = maskChunk && *subRegionMask;
     333             :                     }
     334     7876684 :                     if (anyTrue(maskChunk)) {
     335     7876684 :                         auto chunk = plane(planeBlc, planeTrc);
     336     7876684 :                         if (allTrue(maskChunk)) {
     337     7876659 :                             alg->setData(chunk.begin(), chunk.size());
     338             :                         }
     339             :                         else {
     340          25 :                             alg->setData(chunk.begin(), maskChunk.begin(), chunk.size());
     341             :                         }
     342     7876684 :                         res = alg->getStatistic(_statType);
     343     7876684 :                     }
     344             :                     else {
     345             :                         // no pixels in region on which to do stats, mask grid point
     346           0 :                         writeTo.pixelMask().putAt(False, outPos);
     347             :                     }
     348     7876684 :                     writeTo.putAt(res, outPos);
     349     7876684 :                 }
     350             :             }
     351       42211 :             nCount += nxpts;
     352       42211 :             pm.update(nCount);
     353       42211 :         }
     354         175 :         if (hasLoopAxes) {
     355         174 :             Bool done = False;
     356         174 :             uInt idx = 0;
     357         516 :             while (! done) {
     358         342 :                 auto targetAxis = loopAxes[idx];
     359         342 :                 ++blc[targetAxis];
     360         342 :                 if (blc[targetAxis] == imshape[targetAxis]) {
     361         201 :                     blc[targetAxis] = 0;
     362         201 :                     ++idx;
     363         201 :                     done = idx == loopAxes.size();
     364             :                 }
     365             :                 else {
     366         141 :                     done = True;
     367             :                 }
     368             :             }
     369             :         }
     370         175 :     }
     371          34 : }
     372             : 
     373             : std::shared_ptr<
     374             :     StatisticsAlgorithm<Double, Array<Float>::const_iterator,
     375             :     Array<Bool>::const_iterator>
     376             : >
     377          34 : StatImageCreator::_getStatsAlgorithm(String& algName) const {
     378          34 :     auto ac = _getAlgConf();
     379          34 :     switch(ac.algorithm) {
     380          29 :     case StatisticsData::CLASSICAL:
     381          29 :         algName = "classical";
     382          58 :         return std::shared_ptr<
     383             :             ClassicalStatistics<Double, Array<Float>::const_iterator,
     384             :             Array<Bool>::const_iterator>
     385             :         >(
     386             :             new ClassicalStatistics<
     387             :                 Double, Array<Float>::const_iterator,
     388             :                 Array<Bool>::const_iterator
     389          29 :             >()
     390          29 :         );
     391           5 :     case StatisticsData::CHAUVENETCRITERION:
     392           5 :         algName = "Chauvenet criterion/z-score";
     393          10 :         return std::shared_ptr<ChauvenetCriterionStatistics<
     394             :             Double, Array<Float>::const_iterator,
     395             :             Array<Bool>::const_iterator>
     396             :         >(
     397             :             new ChauvenetCriterionStatistics<
     398             :                 Double, Array<Float>::const_iterator,
     399             :                 Array<Bool>::const_iterator
     400             :             >(
     401             :                 ac.zs, ac.mi
     402           5 :             )
     403           5 :         );
     404           0 :     default:
     405           0 :         ThrowCc("Unsupported statistics algorithm");
     406             :     }
     407             : }
     408             : 
     409          34 : void StatImageCreator::_nominalChunkInfo(
     410             :     std::shared_ptr<Array<Bool>>& templateMask, uInt& xBlcOff, uInt& yBlcOff,
     411             :     uInt& xChunkSize, uInt& yChunkSize, SPCIIF subimage
     412             : ) const {
     413          34 :     Double xPixLength = 0;
     414          34 :     const auto& csys = subimage->coordinates();
     415          34 :     if (_xlen.getUnit().startsWith("pix")) {
     416          34 :         xPixLength = _xlen.getValue();
     417             :     }
     418             :     else {
     419           0 :         const auto& dc = csys.directionCoordinate();
     420           0 :         auto units = dc.worldAxisUnits();
     421           0 :         auto inc = dc.increment();
     422           0 :         Quantity worldPixSize(abs(inc[0]), units[0]);
     423           0 :         xPixLength = _xlen.getValue("rad")/worldPixSize.getValue("rad");
     424           0 :     }
     425          34 :     const auto shape = subimage->shape();
     426          34 :     ThrowIf(
     427             :         xPixLength >= shape[_dirAxes[0]] - 4,
     428             :         "x region length is nearly as large as or larger than the subimage extent of "
     429             :         + String::toString(shape[_dirAxes[0]]) + " pixels. Such a configuration is not supported"
     430             :     );
     431          34 :     ThrowIf(
     432             :         xPixLength <= 0.5,
     433             :         "x region length is only " + String::toString(xPixLength)
     434             :         + " pixels. Such a configuration is not supported"
     435             :     );
     436          34 :     Double yPixLength = xPixLength;
     437          34 :     if (_ylen.getValue() > 0) {
     438          33 :         if (_ylen.getUnit().startsWith("pix")) {
     439          33 :             yPixLength = _ylen.getValue();
     440             :         }
     441             :         else {
     442           0 :             const auto& dc = csys.directionCoordinate();
     443           0 :             auto units = dc.worldAxisUnits();
     444           0 :             auto inc = dc.increment();
     445           0 :             Quantity worldPixSize(abs(inc[1]), units[1]);
     446           0 :             yPixLength = _ylen.getValue("rad")/worldPixSize.getValue("rad");
     447           0 :         }
     448             :     }
     449          34 :     ThrowIf(
     450             :         yPixLength >= shape[_dirAxes[1]] - 4,
     451             :         "y region length is nearly as large as or larger than the subimage extent of "
     452             :         + String::toString(_dirAxes[1]) + " pixels. Such a configuration is not supported"
     453             :     );
     454          34 :     ThrowIf(
     455             :         yPixLength <= 0.5,
     456             :         "y region length is only " + String::toString(yPixLength)
     457             :         + " pixels. Such a configuration is not supported"
     458             :     );
     459          34 :     IPosition templateShape(shape.size(), 1);
     460          34 :     templateShape[_dirAxes[0]] = shape[_dirAxes[0]];
     461          34 :     templateShape[_dirAxes[1]] = shape[_dirAxes[1]];
     462          34 :     TempImage<Float> templateImage(templateShape, csys);
     463             :     // integer division, xq, yq must have integral values
     464          34 :     uInt xcenter = templateShape[_dirAxes[0]]/2;
     465          34 :     uInt ycenter = templateShape[_dirAxes[1]]/2;
     466          34 :     Quantity xq(xcenter, "pix");
     467          34 :     Quantity yq(ycenter, "pix");
     468          34 :     IPosition centerPix(templateShape.size(), 0);
     469          34 :     centerPix[_dirAxes[0]] = xcenter;
     470          34 :     centerPix[_dirAxes[1]] = ycenter;
     471          34 :     auto world = csys.toWorld(centerPix);
     472          34 :     static const Vector<Stokes::StokesTypes> dummyStokes;
     473          34 :     Record reg;
     474          34 :     if (_ylen.getValue() > 0) {
     475             :         AnnCenterBox box(
     476          33 :             xq, yq, _xlen, _ylen, csys, templateShape, dummyStokes
     477          33 :         );
     478          33 :         reg = box.getRegion()->toRecord("");
     479          33 :     }
     480             :     else {
     481             :         AnnCircle circle(
     482           1 :             xq, yq, _xlen, csys, templateShape, dummyStokes
     483           1 :         );
     484           1 :         reg = circle.getRegion()->toRecord("");
     485           1 :     }
     486          34 :     static const AxesSpecifier dummyAxesSpec;
     487             :     auto chunkImage = SubImageFactory<Float>::createSubImageRO(
     488             :         templateImage, reg, "", nullptr, dummyAxesSpec, False
     489          34 :     );
     490          34 :     auto blcOff = chunkImage->coordinates().toPixel(world);
     491          34 :     xBlcOff = rint(blcOff[_dirAxes[0]]);
     492          34 :     yBlcOff = rint(blcOff[_dirAxes[1]]);
     493          34 :     auto chunkShape = chunkImage->shape();
     494          34 :     xChunkSize = chunkShape[_dirAxes[0]];
     495          34 :     yChunkSize = chunkShape[_dirAxes[1]];
     496          34 :     templateMask.reset();
     497          34 :     if (chunkImage->isMasked()) {
     498          34 :         auto mask = chunkImage->getMask();
     499          34 :         if (! allTrue(mask)) {
     500           1 :             templateMask.reset(new Array<Bool>(mask));
     501             :         }
     502          34 :     }
     503          34 : }
     504             : 
     505           6 : void StatImageCreator::_doInterpolation(
     506             :     SPIIF output, TempImage<Float>& store,
     507             :     SPCIIF subImage, uInt nxpts, uInt nypts,
     508             :     Int xstart, Int ystart
     509             : ) const {
     510           6 :     const auto imshape = subImage->shape();
     511           6 :     const auto xshape = imshape[_dirAxes[0]];
     512           6 :     const auto yshape = imshape[_dirAxes[1]];
     513           6 :     const auto ndim = subImage->ndim();
     514          12 :     *_getLog() << LogOrigin(getClass(), __func__)
     515             :         << LogIO::NORMAL << "Interpolate using "
     516           6 :         << _interpName << " algorithm." << LogIO::POST;
     517           6 :     Matrix<Float> result(xshape, yshape);
     518           6 :     Matrix<Bool> resultMask;
     519           6 :     if (_doMask) {
     520           0 :         resultMask.resize(IPosition(2, xshape, yshape));
     521           0 :         resultMask.set(True);
     522             :     }
     523           6 :     Matrix<Float> storage(nxpts, nypts);
     524           6 :     Matrix<Bool> storeMask(nxpts, nypts);
     525           6 :     std::pair<uInt, uInt> start;
     526           6 :     start.first = xstart;
     527           6 :     start.second = ystart;
     528           6 :     if (ndim == 2) {
     529           1 :         store.get(storage);
     530           1 :         store.getMask(storeMask);
     531           1 :         _interpolate(result, resultMask, storage, storeMask, start);
     532           1 :         output->put(result);
     533           1 :         if (_doMask) {
     534           0 :             output->pixelMask().put(resultMask && subImage->pixelMask().get());
     535             :         }
     536             :     }
     537             :     else {
     538             :         // get each direction plane in the storage lattice at each chan/stokes
     539           5 :         IPosition cursorShape(ndim, 1);
     540           5 :         cursorShape[_dirAxes[0]] = nxpts;
     541           5 :         cursorShape[_dirAxes[1]] = nypts;
     542           5 :         auto axisPath = _dirAxes;
     543           5 :         axisPath.append((IPosition::otherAxes(ndim, _dirAxes)));
     544           5 :         LatticeStepper stepper(store.shape(), cursorShape, axisPath);
     545          10 :         Slicer slicer(stepper.position(), stepper.endPosition(), Slicer::endIsLast);
     546           5 :         IPosition myshape(ndim, 1);
     547           5 :         myshape[_dirAxes[0]] = xshape;
     548           5 :         myshape[_dirAxes[1]] = yshape;
     549          12 :         for (stepper.reset(); ! stepper.atEnd(); stepper++) {
     550           7 :             auto pos = stepper.position();
     551           7 :             slicer.setStart(pos);
     552           7 :             slicer.setEnd(stepper.endPosition());
     553           7 :             storage = store.getSlice(slicer, True);
     554           7 :             storeMask = store.getMaskSlice(slicer, True);
     555           7 :             _interpolate(result, resultMask, storage, storeMask, start);
     556           7 :             output->putSlice(result, pos);
     557           7 :             if (_doMask) {
     558           0 :                 output->pixelMask().putSlice(
     559           0 :                     resultMask && subImage->pixelMask().getSlice(pos, myshape, True),
     560             :                     pos
     561             :                 );
     562             :             }
     563           7 :         }
     564           5 :     }
     565           6 : }
     566             : 
     567          34 : void StatImageCreator::setInterpAlgorithm(Interpolate2D::Method alg) {
     568          34 :     switch (alg) {
     569          33 :     case Interpolate2D::CUBIC:
     570          33 :         _interpName = "CUBIC";
     571          33 :         break;
     572           0 :     case Interpolate2D::LANCZOS:
     573           0 :         _interpName = "LANCZOS";
     574           0 :         break;
     575           1 :     case Interpolate2D::LINEAR:
     576           1 :         _interpName = "LINEAR";
     577           1 :         break;
     578           0 :     case Interpolate2D::NEAREST:
     579           0 :         _interpName = "NEAREST";
     580           0 :         break;
     581           0 :     default:
     582           0 :         ThrowCc("Unhandled interpolation method " + String::toString(alg));
     583             :     }
     584          34 :     _interpolater = Interpolate2D(alg);
     585          34 : }
     586             : 
     587          34 : void StatImageCreator::setStatType(StatisticsData::STATS s) {
     588          34 :     _statType = s;
     589          34 :     _statName = StatisticsData::toString(s);
     590          34 :     _statName.upcase();
     591          34 : }
     592             : 
     593          34 : void StatImageCreator::setStatType(const String& s) {
     594          34 :     String m = s;
     595             :     StatisticsData::STATS stat;
     596          34 :     m.downcase();
     597          34 :     if (m.startsWith("i")) {
     598           0 :         stat = StatisticsData::INNER_QUARTILE_RANGE;
     599             :     }
     600          34 :     else if (m.startsWith("max")) {
     601           0 :         stat = StatisticsData::MAX;
     602             :     }
     603          34 :     else if (m.startsWith("mea")) {
     604           0 :         stat = StatisticsData::MEAN;
     605             :     }
     606          34 :     else if (m.startsWith("meda") || m.startsWith("mad") || m.startsWith("x")) {
     607           0 :         stat = StatisticsData::MEDABSDEVMED;
     608             :     }
     609          34 :     else if (m.startsWith("medi")) {
     610           1 :         stat = StatisticsData::MEDIAN;
     611             :     }
     612          33 :     else if (m.startsWith("mi")) {
     613           0 :         stat = StatisticsData::MIN;
     614             :     }
     615          33 :     else if (m.startsWith("n")) {
     616           1 :         stat = StatisticsData::NPTS;
     617             :     }
     618          32 :     else if (m.startsWith("q1")) {
     619           0 :         stat = StatisticsData::FIRST_QUARTILE;
     620             :     }
     621          32 :     else if (m.startsWith("q3")) {
     622           0 :         stat = StatisticsData::THIRD_QUARTILE;
     623             :     }
     624          32 :     else if (m.startsWith("r")) {
     625           0 :         stat = StatisticsData::RMS;
     626             :     }
     627          32 :     else if (m.startsWith("si") || m.startsWith("st")) {
     628          31 :         stat = StatisticsData::STDDEV;
     629             :     }
     630           1 :     else if (m.startsWith("sums")) {
     631           0 :         stat = StatisticsData::SUMSQ;
     632             :     }
     633           1 :     else if (m.startsWith("sum")) {
     634           1 :         stat = StatisticsData::SUM;
     635             :     }
     636           0 :     else if (m.startsWith("v")) {
     637           0 :         stat = StatisticsData::VARIANCE;
     638             :     }
     639             :     else {
     640           0 :         ThrowCc("Statistic " + s + " not supported.");
     641             :     }
     642          34 :     _doProbit = m.startsWith("x");
     643          34 :     setStatType(stat);
     644          34 : }
     645             : 
     646           8 : void StatImageCreator::_interpolate(
     647             :     Matrix<Float>& result, Matrix<Bool>& resultMask,
     648             :     const Matrix<Float>& storage,
     649             :     const Matrix<Bool>& storeMask,
     650             :     const std::pair<uInt, uInt>& start
     651             : ) const {
     652           8 :     auto shape = result.shape();
     653           8 :     auto imax = shape[0];
     654           8 :     auto jmax = shape[1];
     655             :     // x values change fastest in the iterator
     656           8 :     auto iter = result.begin();
     657           8 :     auto miter = resultMask.begin();
     658             :     // xcell and ycell are the current positions within the current
     659             :     // grid cell represented by an integer modulo pointsPerCell
     660           8 :     auto xCell = start.first == 0 ? 0 : _grid.first - start.first;
     661           8 :     auto yCell = start.second == 0 ? 0 : _grid.second - start.second;
     662           8 :     Int xStoreInt = start.first == 0 ? 0 : -1;
     663           8 :     Int yStoreInt = start.second == 0 ? 0 : -1;
     664           8 :     Double xStoreFrac = (Double)xCell/(Double)_grid.first;
     665           8 :     Double yStoreFrac = (Double)yCell/(Double)_grid.second;
     666           8 :     Vector<Double> storeLoc(2);
     667           8 :     storeLoc[0] = xStoreInt + xStoreFrac;
     668           8 :     storeLoc[1] = yStoreInt + yStoreFrac;
     669        1408 :     for (uInt j=0; j<jmax; ++j, ++yCell) {
     670        1400 :         Bool onYGrid = yCell == 0;
     671        1400 :         if (yCell == _grid.second) {
     672         404 :             yCell = 0;
     673         404 :             ++yStoreInt;
     674         404 :             yStoreFrac = 0;
     675         404 :             onYGrid = True;
     676             :         }
     677             :         else {
     678         996 :             yStoreFrac = (Double)yCell/(Double)_grid.second;
     679             :         }
     680        1400 :         storeLoc[1] = yStoreInt + yStoreFrac;
     681        1400 :         xCell = start.first == 0 ? 0 : _grid.first - start.first;
     682        1400 :         xStoreInt = start.first == 0 ? 0 : -1;
     683        1400 :         xStoreFrac = (Double)xCell/(Double)_grid.first;
     684      321400 :         for (uInt i=0; i<imax; ++i, ++xCell) {
     685      320000 :             Bool onXGrid = xCell == 0;
     686      320000 :             if (xCell == _grid.first) {
     687      103100 :                 xCell = 0;
     688      103100 :                 ++xStoreInt;
     689      103100 :                 onXGrid = True;
     690             :             }
     691      320000 :             if (onXGrid && onYGrid) {
     692             :                 // exactly on a grid point, no interpolation needed
     693             :                 // just copy value directly from storage matrix
     694       36149 :                 *iter = storage(xStoreInt, yStoreInt);
     695       36149 :                 if (_doMask) {
     696           0 :                     *miter = storeMask(xStoreInt, yStoreInt);
     697             :                 }
     698             :             }
     699             :             else {
     700      283851 :                 xStoreFrac = (Double)xCell/(Double)_grid.first;
     701      283851 :                 storeLoc[0] = xStoreInt + xStoreFrac;
     702      283851 :                 if (! _interpolater.interp(*iter, storeLoc, storage, storeMask)) {
     703           0 :                     ThrowIf(
     704             :                         ! _doMask,
     705             :                         "Logic error: bad interpolation but there is no mask to set"
     706             :                     );
     707           0 :                     *miter = False;
     708             :                 }
     709             :             }
     710      320000 :             ++iter;
     711      320000 :             if (_doMask) {
     712           0 :                 ++miter;
     713             :             }
     714             :         }
     715             :     }
     716           8 : }
     717             : 
     718             : }