Skip to content

Commit

Permalink
Add more flexible mask propagation to statisticsStack.
Browse files Browse the repository at this point in the history
  • Loading branch information
@TallJimbo
TallJimbo committed Jan 29, 2018
1 parent 2835ff5 commit 7cd263c
Show file tree
Hide file tree
Showing 4 changed files with 165 additions and 18 deletions.
67 changes: 64 additions & 3 deletions include/lsst/afw/math/Stack.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,25 +67,86 @@ void statisticsStack(
std::vector<lsst::afw::image::VariancePixel>(0) ///< vector containing weights
);


/**
* A function to compute some statistics of a stack of Masked Images
*
* @param[in] images MaskedImages to process.
* @param[in] flags Statistics requested.
* @param[in] sctrl Control structure.
* @param[in] wvector Vector of weights.
* @param[in] clipped Mask to set for pixels that were clipped (NOT rejected
* due to masks).
* @param[in] maskMap Vector of pairs of mask pixel values; any pixel
* on an input with any of the bits in .first will result
* in all of the bits in .second being set on the
* corresponding pixel on the output.
*
* If none of the input images are valid for some pixel,
* the afwMath::StatisticsControl::getNoGoodPixelsMask() bit(s) are set.
*
* All the work is done in the function computeMaskedImageStack.
*/
template <typename PixelT>
std::shared_ptr<lsst::afw::image::MaskedImage<PixelT>> statisticsStack(
std::vector<std::shared_ptr<lsst::afw::image::MaskedImage<PixelT>>>& images,
Property flags,
StatisticsControl const& sctrl,
std::vector<lsst::afw::image::VariancePixel> const& wvector,
image::MaskPixel clipped,
std::vector<std::pair<image::MaskPixel, image::MaskPixel>> const & maskMap
);

/**
* A function to compute some statistics of a stack of Masked Images
*
* If none of the input images are valid for some pixel,
* the afwMath::StatisticsControl::getNoGoodPixelsMask() bit(s) are set.
*
* Delegates to the more general version of statisticsStack taking a maskMap.
*/
template <typename PixelT>
std::shared_ptr<lsst::afw::image::MaskedImage<PixelT>> statisticsStack(
std::vector<std::shared_ptr<lsst::afw::image::MaskedImage<PixelT>>>&
images, ///< MaskedImages to process
Property flags, ///< statistics requested
StatisticsControl const& sctrl = StatisticsControl(), ///< control structure
std::vector<lsst::afw::image::VariancePixel> const& wvector =
std::vector<lsst::afw::image::VariancePixel>(0), ///< vector containing weights
image::MaskPixel clipped=0, ///< bitmask to set if any input was clipped
image::MaskPixel clipped=0, ///< bitmask to set if any input was clipped or masked
image::MaskPixel excuse=0 ///< bitmask to excuse from marking as clipped
);
);


/**
* A function to compute some statistics of a stack of Masked Images
*
* @param[out] out Output MaskedImage.
* @param[in] images MaskedImages to process.
* @param[in] flags Statistics requested.
* @param[in] sctrl Control structure.
* @param[in] wvector Vector of weights.
* @param[in] clipped Mask to set for pixels that were clipped (NOT rejected
* due to masks).
* @param[in] maskMap Vector of pairs of mask pixel values; any pixel
* on an input with any of the bits in .first will result
* in all of the bits in .second being set on the
* corresponding pixel on the output.
*
* If none of the input images are valid for some pixel,
* the afwMath::StatisticsControl::getNoGoodPixelsMask() bit(s) are set.
*
* All the work is done in the function computeMaskedImageStack.
*/
template <typename PixelT>
void statisticsStack(lsst::afw::image::MaskedImage<PixelT>& out,
std::vector<std::shared_ptr<lsst::afw::image::MaskedImage<PixelT>>>& images,
Property flags,
StatisticsControl const& sctrl,
std::vector<lsst::afw::image::VariancePixel> const& wvector,
image::MaskPixel clipped,
std::vector<std::pair<image::MaskPixel, image::MaskPixel>> const & maskMap
);

/**
* @ brief compute statistical stack of MaskedImage. Write to output image in-situ
Expand All @@ -98,7 +159,7 @@ void statisticsStack(lsst::afw::image::MaskedImage<PixelT>& out, ///< Output im
StatisticsControl const& sctrl = StatisticsControl(), ///< control structure
std::vector<lsst::afw::image::VariancePixel> const& wvector =
std::vector<lsst::afw::image::VariancePixel>(0), ///< vector containing weights
image::MaskPixel clipped=0, ///< bitmask to set if any input was clipped
image::MaskPixel clipped=0, ///< bitmask to set if any input was clipped or masked
image::MaskPixel excuse=0 ///< bitmask to excuse from marking as clipped
);

Expand Down
18 changes: 18 additions & 0 deletions python/lsst/afw/math/stack.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,15 @@ void declareStatisticsStack(py::module &mod) {
"out"_a, "images"_a, "flags"_a, "sctrl"_a = StatisticsControl(),
"wvector"_a = std::vector<lsst::afw::image::VariancePixel>(0),
"clipped"_a=0, "excuse"_a=0);
mod.def("statisticsStack",
(void (*)(lsst::afw::image::MaskedImage<PixelT> &,
std::vector<std::shared_ptr<lsst::afw::image::MaskedImage<PixelT>>> &, Property,
StatisticsControl const &,
std::vector<lsst::afw::image::VariancePixel> const &,
lsst::afw::image::MaskPixel,
std::vector<std::pair<lsst::afw::image::MaskPixel, lsst::afw::image::MaskPixel>> const &
))statisticsStack<PixelT>,
"out"_a, "images"_a, "flags"_a, "sctrl"_a, "wvector"_a, "clipped"_a, "maskMap"_a);
mod.def("statisticsStack",
(std::shared_ptr<lsst::afw::image::Image<PixelT>>(*)(
std::vector<std::shared_ptr<lsst::afw::image::Image<PixelT>>> &, Property,
Expand All @@ -74,6 +83,15 @@ void declareStatisticsStack(py::module &mod) {
"images"_a, "flags"_a, "sctrl"_a = StatisticsControl(),
"wvector"_a = std::vector<lsst::afw::image::VariancePixel>(0),
"clipped"_a=0, "excuse"_a=0);
mod.def("statisticsStack",
(std::shared_ptr<lsst::afw::image::MaskedImage<PixelT>>(*)(
std::vector<std::shared_ptr<lsst::afw::image::MaskedImage<PixelT>>> &, Property,
StatisticsControl const &,
std::vector<lsst::afw::image::VariancePixel> const &,
lsst::afw::image::MaskPixel,
std::vector<std::pair<lsst::afw::image::MaskPixel, lsst::afw::image::MaskPixel>> const &
))statisticsStack<PixelT>,
"images"_a, "flags"_a, "sctrl"_a, "wvector"_a, "clipped"_a, "maskMap"_a);
mod.def("statisticsStack",
(std::shared_ptr<std::vector<PixelT>>(*)(
std::vector<std::shared_ptr<std::vector<PixelT>>> &, Property, StatisticsControl const &,
Expand Down
87 changes: 72 additions & 15 deletions src/math/Stack.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -102,6 +102,8 @@ void checkImageSizes(ImageT const &out, std::vector<std::shared_ptr<ImageT>> con
* stack MaskedImages
*
* ************************************************************************** */

//@{
/**
* @internal A function to handle MaskedImage stacking
*
Expand All @@ -115,8 +117,10 @@ void computeMaskedImageStack(image::MaskedImage<PixelT> &imgStack,
std::vector<std::shared_ptr<image::MaskedImage<PixelT>>> const &images,
Property flags, StatisticsControl const &sctrl,
image::MaskPixel const clipped,
image::MaskPixel const excuse,
std::vector<std::pair<image::MaskPixel, image::MaskPixel>> const & maskMap,
WeightVector const &wvector = WeightVector()) {


// get a list of row_begin iterators
typedef typename image::MaskedImage<PixelT>::x_iterator x_iterator;
std::vector<x_iterator> rows;
Expand Down Expand Up @@ -178,33 +182,41 @@ void computeMaskedImageStack(image::MaskedImage<PixelT> &imgStack,
* getting a variance of NaN when you only have one input
*/
}
// Check to see if any pixels were rejected
// Check to see if any pixels were rejected due to clipping
if (stat.getValue(NCLIPPED) > 0) {
msk |= clipped;
} else if (stat.getValue(NMASKED) > 0) {
bool maskIt = false;
if (excuse != 0) {
// Might the masking be excused?
image::MaskPixel const maskVal = sctrl.getAndMask() & ~excuse;
}
// Check to see if any pixels were rejected to be masking, and apply
// any requested masks to the result.
if (stat.getValue(NMASKED) > 0) {
for (auto const & pair : maskMap) {
for (auto pp = pixelSet.begin(); pp != pixelSet.end(); ++pp) {
if ((*pp).mask() & maskVal) {
maskIt = true;
if ((*pp).mask() & pair.first) {
msk |= pair.second;
break;
}
}
} else {
maskIt = true;
}

if (maskIt) {
msk |= clipped;
}
}

*ptr = typename image::MaskedImage<PixelT>::Pixel(stat.getValue(flags), msk, variance);
}
}
}
template <typename PixelT, bool isWeighted, bool useVariance>
void computeMaskedImageStack(image::MaskedImage<PixelT> &imgStack,
std::vector<std::shared_ptr<image::MaskedImage<PixelT>>> const &images,
Property flags, StatisticsControl const &sctrl,
image::MaskPixel const clipped,
image::MaskPixel const excuse,
WeightVector const &wvector = WeightVector()) {
std::vector<std::pair<image::MaskPixel, image::MaskPixel>> maskMap;
maskMap.push_back(std::make_pair(sctrl.getAndMask() & ~excuse, clipped));
computeMaskedImageStack<PixelT, isWeighted, useVariance>(
imgStack, images, flags, sctrl, clipped, maskMap, wvector
);
}
//@}

} // end anonymous namespace

Expand All @@ -221,6 +233,21 @@ std::shared_ptr<image::MaskedImage<PixelT>> statisticsStack(
statisticsStack(*out, images, flags, sctrl, wvector, clipped, excuse);
return out;
}

template <typename PixelT>
std::shared_ptr<image::MaskedImage<PixelT>> statisticsStack(
std::vector<std::shared_ptr<image::MaskedImage<PixelT>>> &images, Property flags,
StatisticsControl const &sctrl, WeightVector const &wvector, image::MaskPixel clipped,
std::vector<std::pair<image::MaskPixel, image::MaskPixel>> const & maskMap) {
if (images.size() == 0) {
throw LSST_EXCEPT(pexExcept::LengthError, "Please specify at least one image to stack");
}
std::shared_ptr<image::MaskedImage<PixelT>> out(
new image::MaskedImage<PixelT>(images[0]->getDimensions()));
statisticsStack(*out, images, flags, sctrl, wvector, clipped, maskMap);
return out;
}

template <typename PixelT>
void statisticsStack(image::MaskedImage<PixelT> &out,
std::vector<std::shared_ptr<image::MaskedImage<PixelT>>> &images, Property flags,
Expand All @@ -243,6 +270,28 @@ void statisticsStack(image::MaskedImage<PixelT> &out,
}
}

template <typename PixelT>
void statisticsStack(image::MaskedImage<PixelT> &out,
std::vector<std::shared_ptr<image::MaskedImage<PixelT>>> &images, Property flags,
StatisticsControl const &sctrl, WeightVector const &wvector, image::MaskPixel clipped,
std::vector<std::pair<image::MaskPixel, image::MaskPixel>> const & maskMap) {
checkObjectsAndWeights(images, wvector);
checkOnlyOneFlag(flags);
checkImageSizes(out, images);

if (sctrl.getWeighted()) {
if (wvector.empty()) {
return computeMaskedImageStack<PixelT, true, true>(out, images, flags, sctrl,
clipped, maskMap); // use variance
} else {
return computeMaskedImageStack<PixelT, true, false>(out, images, flags, sctrl, clipped, maskMap,
wvector); // use wvector
}
} else {
return computeMaskedImageStack<PixelT, false, false>(out, images, flags, sctrl, clipped, maskMap);
}
}

namespace {
/* ************************************************************************** *
*
Expand Down Expand Up @@ -487,10 +536,18 @@ std::shared_ptr<image::MaskedImage<PixelT>> statisticsStack(image::MaskedImage<P
std::vector<std::shared_ptr<image::MaskedImage<TYPE>>> & images, Property flags, \
StatisticsControl const &sctrl, WeightVector const &wvector, image::MaskPixel, \
image::MaskPixel); \
template std::shared_ptr<image::MaskedImage<TYPE>> statisticsStack<TYPE>( \
std::vector<std::shared_ptr<image::MaskedImage<TYPE>>> & images, Property flags, \
StatisticsControl const &sctrl, WeightVector const &wvector, image::MaskPixel, \
std::vector<std::pair<image::MaskPixel, image::MaskPixel>> const &); \
template void statisticsStack<TYPE>( \
image::MaskedImage<TYPE> & out, std::vector<std::shared_ptr<image::MaskedImage<TYPE>>> & images, \
Property flags, StatisticsControl const &sctrl, WeightVector const &wvector, image::MaskPixel, \
image::MaskPixel); \
template void statisticsStack<TYPE>( \
image::MaskedImage<TYPE> & out, std::vector<std::shared_ptr<image::MaskedImage<TYPE>>> & images, \
Property flags, StatisticsControl const &sctrl, WeightVector const &wvector, image::MaskPixel, \
std::vector<std::pair<image::MaskPixel, image::MaskPixel>> const &); \
template std::shared_ptr<std::vector<TYPE>> statisticsStack<TYPE>( \
std::vector<std::shared_ptr<std::vector<TYPE>>> & vectors, Property flags, \
StatisticsControl const &sctrl, WeightVector const &wvector); \
Expand Down
11 changes: 11 additions & 0 deletions tests/test_stacker.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -398,6 +398,17 @@ def testClipped(self):
self.assertFloatsAlmostEqual(stack.getImage().getArray(), 0.0, atol=0.0)
self.assertEqual(stack.getMask().get(1, 1), 0)

# Map that mask value to a different one.
rejected = 1 << afwImage.Mask().addMaskPlane("REJECTED")
maskMap = [(maskVal, rejected)]
images[0].getMask().set(1, 1, 0) # only want to clip, not mask, this one
images[1].getMask().set(1, 2, maskVal) # only want to mask, not clip, this one
stack = afwMath.statisticsStack(images, afwMath.MEANCLIP, statsCtrl, wvector=[], clipped=clipped,
maskMap=maskMap)
self.assertFloatsAlmostEqual(stack.getImage().getArray(), 0.0, atol=0.0)
self.assertEqual(stack.getMask().get(1, 1), clipped)
self.assertEqual(stack.getMask().get(1, 2), rejected)

#################################################################
# Test suite boiler plate
#################################################################
Expand Down

0 comments on commit 7cd263c

Please sign in to comment.