Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

DM-13396: Add more flexible mask propagation to statisticsStack. #314

Merged
merged 1 commit into from
Jan 30, 2018
Merged

DM-13396: Add more flexible mask propagation to statisticsStack. #314

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
67 changes: 64 additions & 3 deletions include/lsst/afw/math/Stack.h
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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 by masking, and apply
// any associated 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
Open in desktop
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