casa / casa5 / 6541319b7b5

Commits

Ville Suoranta authored 6541319b7b5 Merge 13 Jun 2017

Merge pull request #99 in CASA/casa from bugfix/CAS-10109 to master

* commit 'f6e3baa42a60c1fdf2ef973c6e61c872c6742f1a':
  added docs regarding masking of output image
  mask output pixels where interpolation was unsuccessful
  implement copying of input mask to output

code/imageanalysis/ImageAnalysis/StatImageCreator.cc

Modified

-64
-65
+void StatImageCreator::setGridSpacing(uInt x, uInt y) {
-66
+    _grid.first = x;
-67
+    _grid.second = y;
-68
+}
-69
-70
+SPIIF StatImageCreator::compute() {
-71
+    static const AxesSpecifier dummyAxesSpec;
-72
+    *_getLog() << LogOrigin(getClass(), __func__);
-73
+    auto mylog = _getLog().get();
--
+    auto subImage = SubImageFactory<Float>::createSubImageRO(
-+
+    auto subImageRO = SubImageFactory<Float>::createSubImageRO(
-75
+        *_getImage(), *_getRegion(), _getMask(),
-76
+        mylog, dummyAxesSpec, _getStretch()
-77
+    );
--
+    _doMask = ! ImageMask::isAllMaskTrue(*subImage);
--
+    const auto imshape = subImage->shape();
-+
+    auto subImageRW = SubImageFactory<Float>::createImage(
-+
+        *_getImage(), "", *_getRegion(), _getMask(),
-+
+        dummyAxesSpec, False, False, _getStretch()
-+
+    );
-+
+    _doMask = ! ImageMask::isAllMaskTrue(*subImageRO);
-+
+    const auto imshape = subImageRO->shape();
-84
+    const auto xshape = imshape[_dirAxes[0]];
-85
+    const auto yshape = imshape[_dirAxes[1]];
--
+    const auto& csys = subImage->coordinates();
-+
+    const auto& csys = subImageRO->coordinates();
-87
+    auto anchorPixel = csys.toPixel(_anchor);
--
+    TempImage<Float> output(imshape, csys);
--
+    output.set(0);
--
+    if (_doMask) {
--
+        output.attachMask(ArrayLattice<Bool>(imshape));
--
+        output.pixelMask().set(True);
--
+    }
-88
+    Int xanchor = rint(anchorPixel[_dirAxes[0]]);
-89
+    Int yanchor = rint(anchorPixel[_dirAxes[1]]);
-90
+    // ensure xanchor and yanchor are positive
-91
+    if (xanchor < 0) {
-92
+        // ugh, mod of a negative number in C++ doesn't do what I want it to
-93
+        // integer division
-94
+        xanchor += (abs(xanchor)/_grid.first + 1)*_grid.first;
-95
+    }
-96
+    if (yanchor < 0) {
-97
+        yanchor += (abs(yanchor)/_grid.second + 1)*_grid.second;
-115
+        ++nxpts;
-116
+    }
-117
+    auto nypts = (yshape - ystart)/ _grid.second;
-118
+    if ((yshape - ystart) % _grid.second != 0) {
-119
+        ++nypts;
-120
+    }
-121
+    IPosition storeShape = imshape;
-122
+    storeShape[_dirAxes[0]] = nxpts;
-123
+    storeShape[_dirAxes[1]] = nypts;
-124
+    auto interpolate = _grid.first > 1 || _grid.second > 1;
--
+    TempImage<Float>* writeTo = &output;
-+
+    TempImage<Float>* writeTo = dynamic_cast<TempImage<Float> *>(subImageRW.get());
-126
+    std::unique_ptr<TempImage<Float>> store;
-127
+    if (interpolate) {
-128
+        store.reset(new TempImage<Float>(storeShape, csys));
-129
+        store->set(0.0);
-130
+        if (_doMask) {
-131
+            store->attachMask(ArrayLattice<Bool>(storeShape));
-132
+            store->pixelMask().set(True);
-133
+        }
-134
+        writeTo = store.get();
-135
+    }
--
+    _computeStat(*writeTo, subImage, nxpts, nypts, xstart, ystart);
-+
+    _computeStat(*writeTo, subImageRO, nxpts, nypts, xstart, ystart);
-137
+    if (_doPhi) {
-138
+        writeTo->copyData((LatticeExpr<Float>)(*writeTo * PHI));
-139
+    }
--
-140
+    if (interpolate) {
-141
+        _doInterpolation(
--
+            output, *store,
--
+            subImage,  nxpts,  nypts,
--
+            xstart,  ystart
-+
+            subImageRW, *store, subImageRO,
-+
+            nxpts, nypts, xstart,  ystart
-144
+        );
-145
+    }
--
+    auto res = _prepareOutputImage(output);
-+
+    auto res = _prepareOutputImage(*subImageRW);
-147
+    return res;
-148
+}
-149
-150
+void StatImageCreator::_computeStat(
-151
+    TempImage<Float>& writeTo, SPCIIF subImage,
-152
+    uInt nxpts, uInt nypts, Int xstart, Int ystart
-153
+) {
-154
+    SHARED_PTR<Array<Bool>> regionMask;
-155
+    uInt xBlcOff = 0;
-156
+    uInt yBlcOff = 0;
-175
+    chunkShape[_dirAxes[0]] = xChunkSize;
-176
+    chunkShape[_dirAxes[1]] = yChunkSize;
-177
+    auto xsize = imshape[_dirAxes[0]];
-178
+    auto ysize = imshape[_dirAxes[1]];
-179
+    IPosition planeShape(imshape.size(), 1);
-180
+    planeShape[_dirAxes[0]] = xsize;
-181
+    planeShape[_dirAxes[1]] = ysize;
-182
+    RO_MaskedLatticeIterator<Float> lattIter (
-183
+        *subImage, planeShape, True
-184
+    );
--
+    // Vector rather than IPosition because blc values can be negative
--
+    Vector<Int> blc(ndim, 0);
--
+    blc[_dirAxes[0]] = xstart - xBlcOff;
--
+    blc[_dirAxes[1]] = ystart - yBlcOff;
-185
+    String algName;
-186
+    auto alg = _getStatsAlgorithm(algName);
-187
+    std::set<StatisticsData::STATS> statToCompute;
-188
+    statToCompute.insert(_statType);
-189
+    alg->setStatsToCalculate(statToCompute);
-190
+    auto ngrid = nxpts*nypts;
-191
+    auto nPts = ngrid;
-192
+    IPosition loopAxes;
-193
+    for (uInt i=0; i<ndim; ++i) {
-194
+        if (i != _dirAxes[0] && i != _dirAxes[1]) {
-195
+            loopAxes.append(IPosition(1, i));
-196
+            nPts *= imshape[i];
-197
+        }
-198
+    }
--
+    auto hasLoopAxes = loopAxes.size() > 0;
--
+    ProgressMeter pm(0, nPts, "Processing stats at grid points");
-199
+    auto doCircle = _ylen.getValue() <= 0;
-200
+    *_getLog() << LogOrigin(getClass(), __func__) << LogIO::NORMAL
-201
+        << "Using ";
-202
+    if (doCircle) {
-203
+        *_getLog() << "circular region of radius " << _xlen;
-204
+    }
-205
+    else {
-206
+        *_getLog() << "rectangular region of specified dimensions " << _xlen
-207
+            << " x " << _ylen;
-208
+    }
-209
+    *_getLog() << " (because of centering and "
-210
+            << "rounding to use whole pixels, actual dimensions of bounding box are "
-211
+            << xChunkSize << " pix x " << yChunkSize << " pix";
-212
+    if (doCircle) {
-213
+        *_getLog() << " and there are " << ntrue(*regionMask)
-214
+            << " good pixels in the circle that are being used";
-215
+    }
-216
+    *_getLog() << ") to choose pixels for computing " << _statName
-217
+        << " using the " << algName << " algorithm around each of "
-218
+        << ngrid << " grid points in " << (nPts/ngrid) << " planes." << LogIO::POST;
-+
+    auto imageChunkShape = chunkShape;
-+
+    _doStatsLoop(
-+
+        writeTo, lattIter, nxpts, nypts, xstart, ystart,
-+
+        xBlcOff, yBlcOff, xChunkSize, yChunkSize, imshape,
-+
+        chunkShape, regionMask, alg, regMaskCopy, loopAxes, nPts
-+
+    );
-+
+}
-+
-+
+void StatImageCreator::_doStatsLoop(
-+
+    TempImage<Float>& writeTo, RO_MaskedLatticeIterator<Float>& lattIter,
-+
+    uInt nxpts, uInt nypts, Int xstart, Int ystart, uInt xBlcOff,
-+
+    uInt yBlcOff, uInt xChunkSize, uInt yChunkSize, const IPosition& imshape,
-+
+    const IPosition& chunkShape, SHARED_PTR<Array<Bool>> regionMask,
-+
+    SHARED_PTR<
-+
+        StatisticsAlgorithm<
-+
+            Double, Array<Float>::const_iterator, Array<Bool>::const_iterator,
-+
+            Array<Float>::const_iterator
-+
+        >
-+
+    >& alg, const Array<Bool>& regMaskCopy, const IPosition& loopAxes, uInt nPts
-+
+) {
-+
+    auto ximshape = imshape[_dirAxes[0]];
-+
+    auto yimshape = imshape[_dirAxes[1]];
-+
+    auto ndim = imshape.size();
-242
+    IPosition planeBlc(ndim, 0);
-243
+    auto& xPlaneBlc = planeBlc[_dirAxes[0]];
-244
+    auto& yPlaneBlc = planeBlc[_dirAxes[1]];
--
+    auto imageChunkShape = chunkShape;
--
+    // pixels at the TRC are included in the statistics
-245
+    auto planeTrc = planeBlc + chunkShape - 1;
-246
+    auto& xPlaneTrc = planeTrc[_dirAxes[0]];
-247
+    auto& yPlaneTrc = planeTrc[_dirAxes[1]];
-248
+    uInt nCount = 0;
--
+    auto ximshape = imshape[_dirAxes[0]];
--
+    auto yimshape = imshape[_dirAxes[1]];
-+
+    auto hasLoopAxes = ! loopAxes.empty();
-+
+    ProgressMeter pm(0, nPts, "Processing stats at grid points");
-+
+    // Vector rather than IPosition because blc values can be negative
-+
+    Vector<Int> blc(ndim, 0);
-+
+    blc[_dirAxes[0]] = xstart - xBlcOff;
-+
+    blc[_dirAxes[1]] = ystart - yBlcOff;
-255
+    for (lattIter.atStart(); ! lattIter.atEnd(); ++lattIter) {
-256
+        auto plane = lattIter.cursor();
-257
+        auto lattMask = lattIter.getMask();
-+
+        auto checkGrid = ! allTrue(lattMask);
-259
+        auto outPos = lattIter.position();
-260
+        auto& xOutPos = outPos[_dirAxes[0]];
-261
+        auto& yOutPos = outPos[_dirAxes[1]];
-+
+        // inPos is the position of the current grid point
-+
+        // This is the position in the current chunk, not the entire lattice
-+
+        auto inPos = plane.shape() - 1;
-+
+        inPos[_dirAxes[0]] = xstart;
-+
+        inPos[_dirAxes[1]] = ystart;
-+
+        auto& xInPos = inPos[_dirAxes[0]];
-+
+        auto& yInPos = inPos[_dirAxes[1]];
-269
+        Int yblc = ystart - yBlcOff;
--
+        for (uInt yCount=0; yCount<nypts; ++yCount, yblc+=_grid.second) {
-+
+        for (uInt yCount=0; yCount<nypts; ++yCount, yblc+=_grid.second, yInPos+=_grid.second) {
-271
+            yOutPos = yCount;
-272
+            yPlaneBlc = max(0, yblc);
-273
+            yPlaneTrc = min(
-274
+                yblc + (Int)yChunkSize - 1, (Int)yimshape - 1
-275
+            );
-276
+            IPosition regMaskStart(ndim, 0);
-277
+            auto& xRegMaskStart = regMaskStart[_dirAxes[0]];
-278
+            auto& yRegMaskStart = regMaskStart[_dirAxes[1]];
-279
+            auto regMaskLength = regMaskStart;
-280
+            auto& xRegMaskLength = regMaskLength[_dirAxes[0]];
-285
+                if (yblc < 0) {
-286
+                    yRegMaskStart = -yblc;
-287
+                    yRegMaskLength += yblc;
-288
+                    yDoMaskSlice = True;
-289
+                }
-290
+                else if (yblc + yChunkSize > yimshape) {
-291
+                    yRegMaskLength = yimshape - yblc;
-292
+                    yDoMaskSlice = True;
-293
+                }
-294
+            }
-+
+            xInPos = xstart;
-296
+            Int xblc = xstart - xBlcOff;
--
+            for (uInt xCount=0; xCount<nxpts; ++xCount, xblc+=_grid.first) {
-+
+            for (uInt xCount=0; xCount<nxpts; ++xCount, xblc+=_grid.first, xInPos+=_grid.first) {
-+
+                Float res = 0;
-299
+                xOutPos = xCount;
--
+                xPlaneBlc = max(0, xblc);
--
+                xPlaneTrc = min(
--
+                    xblc + (Int)xChunkSize - 1, (Int)ximshape - 1
--
+                );
--
+                SHARED_PTR<Array<Bool>> subRegionMask;
--
+                if (regionMask) {
--
+                    auto doMaskSlice = yDoMaskSlice;
--
+                    xRegMaskStart = 0;
--
+                    if (xblc < 0) {
--
+                        xRegMaskStart = -xblc;
--
+                        xRegMaskLength = regionMask->shape()[_dirAxes[0]] + xblc;
--
+                        doMaskSlice = True;
--
+                    }
--
+                    else if (xblc + xChunkSize > ximshape) {
--
+                        regMaskLength[_dirAxes[0]] = ximshape - xblc;
--
+                        doMaskSlice = True;
--
+                    }
--
+                    else {
--
+                        xRegMaskLength = xChunkSize;
--
+                    }
--
+                    if (doMaskSlice) {
--
+                        Slicer sl(regMaskStart, regMaskLength);
--
+                        subRegionMask.reset(new Array<Bool>(regMaskCopy(sl)));
-+
+                if (checkGrid && ! lattMask(inPos)) {
-+
+                    // grid point is masked, no computation should be done
-+
+                    writeTo.putAt(res, outPos);
-+
+                    writeTo.pixelMask().putAt(False, outPos);
-+
+                }
-+
+                else {
-+
+                    xPlaneBlc = max(0, xblc);
-+
+                    xPlaneTrc = min(
-+
+                        xblc + (Int)xChunkSize - 1, (Int)ximshape - 1
-+
+                    );
-+
+                    SHARED_PTR<Array<Bool>> subRegionMask;
-+
+                    if (regionMask) {
-+
+                        auto doMaskSlice = yDoMaskSlice;
-+
+                        xRegMaskStart = 0;
-+
+                        if (xblc < 0) {
-+
+                            xRegMaskStart = -xblc;
-+
+                            xRegMaskLength = regionMask->shape()[_dirAxes[0]] + xblc;
-+
+                            doMaskSlice = True;
-+
+                        }
-+
+                        else if (xblc + xChunkSize > ximshape) {
-+
+                            regMaskLength[_dirAxes[0]] = ximshape - xblc;
-+
+                            doMaskSlice = True;
-+
+                        }
-+
+                        else {
-+
+                            xRegMaskLength = xChunkSize;
-+
+                        }
-+
+                        if (doMaskSlice) {
-+
+                            Slicer sl(regMaskStart, regMaskLength);
-+
+                            subRegionMask.reset(new Array<Bool>(regMaskCopy(sl)));
-+
+                        }
-+
+                        else {
-+
+                            subRegionMask.reset(new Array<Bool>(regMaskCopy));
-+
+                        }
-333
+                    }
--
+                    else {
--
+                        subRegionMask.reset(new Array<Bool>(regMaskCopy));
-+
+                    auto maskChunk = lattMask(planeBlc, planeTrc).copy();
-+
+                    if (subRegionMask) {
-+
+                        maskChunk = maskChunk && *subRegionMask;
-337
+                    }
--
+                }
--
+                auto maskChunk = lattMask(planeBlc, planeTrc).copy();
--
+                if (subRegionMask) {
--
+                    maskChunk = maskChunk && *subRegionMask;
--
+                }
--
+                Float res = 0;
--
+                if (anyTrue(maskChunk)) {
--
+                    auto chunk = plane(planeBlc, planeTrc);
--
+                    if (allTrue(maskChunk)) {
--
+                        alg->setData(chunk.begin(), chunk.size());
-+
+                    if (anyTrue(maskChunk)) {
-+
+                        auto chunk = plane(planeBlc, planeTrc);
-+
+                        if (allTrue(maskChunk)) {
-+
+                            alg->setData(chunk.begin(), chunk.size());
-+
+                        }
-+
+                        else {
-+
+                            alg->setData(chunk.begin(), maskChunk.begin(), chunk.size());
-+
+                        }
-+
+                        res = alg->getStatistic(_statType);
-347
+                    }
-348
+                    else {
--
+                        alg->setData(chunk.begin(), maskChunk.begin(), chunk.size());
-+
+                        // no pixels in region on which to do stats, mask grid point
-+
+                        writeTo.pixelMask().putAt(False, outPos);
-351
+                    }
--
+                    res = alg->getStatistic(_statType);
-+
+                    writeTo.putAt(res, outPos);
-353
+                }
--
+                else {
--
+                    writeTo.pixelMask().putAt(False, outPos);
--
+                }
--
+                writeTo.putAt(res, outPos);
-354
+            }
-355
+            nCount += nxpts;
-356
+            pm.update(nCount);
-357
+        }
-358
+        if (hasLoopAxes) {
-359
+            Bool done = False;
-360
+            uInt idx = 0;
-361
+            while (! done) {
-362
+                auto targetAxis = loopAxes[idx];
-363
+                ++blc[targetAxis];
-500
+    templateMask.reset();
-501
+    if (chunkImage->isMasked()) {
-502
+        auto mask = chunkImage->getMask();
-503
+        if (! allTrue(mask)) {
-504
+            templateMask.reset(new Array<Bool>(mask));
-505
+        }
-506
+    }
-507
+}
-508
-509
+void StatImageCreator::_doInterpolation(
--
+    TempImage<Float>& output, TempImage<Float>& store,
-+
+    SPIIF output, TempImage<Float>& store,
-511
+    SPCIIF subImage, uInt nxpts, uInt nypts,
-512
+    Int xstart, Int ystart
-513
+) const {
-514
+    const auto imshape = subImage->shape();
-515
+    const auto xshape = imshape[_dirAxes[0]];
-516
+    const auto yshape = imshape[_dirAxes[1]];
-517
+    const auto ndim = subImage->ndim();
-518
+    *_getLog() << LogOrigin(getClass(), __func__)
-519
+        << LogIO::NORMAL << "Interpolate using "
-520
+        << _interpName << " algorithm." << LogIO::POST;
-526
+    }
-527
+    Matrix<Float> storage(nxpts, nypts);
-528
+    Matrix<Bool> storeMask(nxpts, nypts);
-529
+    std::pair<uInt, uInt> start;
-530
+    start.first = xstart;
-531
+    start.second = ystart;
-532
+    if (ndim == 2) {
-533
+        store.get(storage);
-534
+        store.getMask(storeMask);
-535
+        _interpolate(result, resultMask, storage, storeMask, start);
--
+        output.put(result);
-+
+        output->put(result);
-537
+        if (_doMask) {
--
+            output.pixelMask().put(resultMask);
-+
+            output->pixelMask().put(resultMask && subImage->pixelMask().get());
-539
+        }
-540
+    }
-541
+    else {
-542
+        // get each direction plane in the storage lattice at each chan/stokes
-543
+        IPosition cursorShape(ndim, 1);
-544
+        cursorShape[_dirAxes[0]] = nxpts;
-545
+        cursorShape[_dirAxes[1]] = nypts;
-546
+        auto axisPath = _dirAxes;
--
+        axisPath.append((IPosition::otherAxes(ndim,_dirAxes)));
-+
+        axisPath.append((IPosition::otherAxes(ndim, _dirAxes)));
-548
+        LatticeStepper stepper(store.shape(), cursorShape, axisPath);
-549
+        Slicer slicer(stepper.position(), stepper.endPosition(), Slicer::endIsLast);
-+
+        IPosition myshape(ndim, 1);
-+
+        myshape[_dirAxes[0]] = xshape;
-+
+        myshape[_dirAxes[1]] = yshape;
-553
+        for (stepper.reset(); ! stepper.atEnd(); stepper++) {
-554
+            auto pos = stepper.position();
-555
+            slicer.setStart(pos);
-556
+            slicer.setEnd(stepper.endPosition());
-557
+            storage = store.getSlice(slicer, True);
-558
+            storeMask = store.getMaskSlice(slicer, True);
-559
+            _interpolate(result, resultMask, storage, storeMask, start);
--
+            output.putSlice(result, pos);
-+
+            output->putSlice(result, pos);
-561
+            if (_doMask) {
--
+                output.pixelMask().putSlice(resultMask, pos);
-+
+                output->pixelMask().putSlice(
-+
+                    resultMask && subImage->pixelMask().getSlice(pos, myshape, True),
-+
+                    pos
-+
+                );
-566
+            }
-567
+        }
-568
+    }
-569
+}
-570
-571
+void StatImageCreator::setInterpAlgorithm(Interpolate2D::Method alg) {
-572
+    switch (alg) {
-573
+    case Interpolate2D::CUBIC:
-574
+        _interpName = "CUBIC";
-575
+        return;
-696
+                // exactly on a grid point, no interpolation needed
-697
+                // just copy value directly from storage matrix
-698
+                *iter = storage(xStoreInt, yStoreInt);
-699
+                if (_doMask) {
-700
+                    *miter = storeMask(xStoreInt, yStoreInt);
-701
+                }
-702
+            }
-703
+            else {
-704
+                xStoreFrac = (Double)xCell/(Double)_grid.first;
-705
+                storeLoc[0] = xStoreInt + xStoreFrac;
--
+                _interpolater.interp(*iter, storeLoc, storage, storeMask);
-+
+                if (! _interpolater.interp(*iter, storeLoc, storage, storeMask)) {
-+
+                    ThrowIf(
-+
+                        ! _doMask,
-+
+                        "Logic error: bad interpolation but there is no mask to set"
-+
+                    );
-+
+                    *miter = False;
-+
+                }
-713
+            }
-714
+            ++iter;
-715
+            if (_doMask) {
-716
+                ++miter;
-717
+            }
-718
+        }
-719
+    }
-720
+}
-721
-722
+}

casa5

Commits

Jira Issues

Add shortcut