saturated.cc

Go to the documentation of this file.

/*
 * Handle saturated pixels when making colour images
 */
#include <cmath>
#include "boost/format.hpp"
#include "lsst/afw/detection.h"
#include "lsst/afw/image/MaskedImage.h"
#include "rgb/Rgb.h"

namespace lsst {
namespace afw {
namespace display {

namespace {
template <typename ImageT>
class SetPixels {
public:
    explicit SetPixels() : _value(0) {}

    void setValue(float value) { _value = value; }

    void operator()(lsst::geom::Point2I const& point, typename ImageT::Pixel& arrayInput) { arrayInput = _value; }

private:
    float _value;
};
}

template <typename ImageT>
void replaceSaturatedPixels(ImageT& rim,      // R image (e.g. i)
                            ImageT& gim,      // G image (e.g. r)
                            ImageT& bim,      // B image (e.g. g)
                            int borderWidth,  // width of border used to estimate colour of saturated regions
                            float saturatedPixelValue  // the brightness of a saturated pixel, once fixed
                            ) {
    int const width = rim.getWidth(), height = rim.getHeight();
    int const x0 = rim.getX0(), y0 = rim.getY0();

    if (width != gim.getWidth() || height != gim.getHeight() || x0 != gim.getX0() || y0 != gim.getY0()) {
        throw LSST_EXCEPT(
                pex::exceptions::InvalidParameterError,
                str(boost::format("R image has different size/origin from G image "
                                  "(%dx%d+%d+%d v. %dx%d+%d+%d") %
                    width % height % x0 % y0 % gim.getWidth() % gim.getHeight() % gim.getX0() % gim.getY0()));
    }
    if (width != bim.getWidth() || height != bim.getHeight() || x0 != bim.getX0() || y0 != bim.getY0()) {
        throw LSST_EXCEPT(
                pex::exceptions::InvalidParameterError,
                str(boost::format("R image has different size/origin from B image "
                                  "(%dx%d+%d+%d v. %dx%d+%d+%d") %
                    width % height % x0 % y0 % bim.getWidth() % bim.getHeight() % bim.getX0() % bim.getY0()));
    }

    bool const useMaxPixel = !std::isfinite(saturatedPixelValue);

    SetPixels<typename ImageT::Image> setR, setG, setB;  // functors used to set pixel values

    // Find all the saturated pixels in any of the three image
    int const npixMin = 1;  // minimum number of pixels in an object
    afw::image::MaskPixel const SAT = rim.getMask()->getPlaneBitMask("SAT");
    detection::Threshold const satThresh(SAT, detection::Threshold::BITMASK);

    detection::FootprintSet sat(*rim.getMask(), satThresh, npixMin);
    sat.merge(detection::FootprintSet(*gim.getMask(), satThresh, npixMin));
    sat.merge(detection::FootprintSet(*bim.getMask(), satThresh, npixMin));
    // go through the list of saturated regions, determining the mean colour of the surrounding pixels
    typedef detection::FootprintSet::FootprintList FootprintList;
    std::shared_ptr<FootprintList> feet = sat.getFootprints();
    for (FootprintList::const_iterator ptr = feet->begin(), end = feet->end(); ptr != end; ++ptr) {
        std::shared_ptr<detection::Footprint> const foot = *ptr;
        auto const bigFoot = std::make_shared<detection::Footprint>(foot->getSpans()->dilated(borderWidth),
                                                                    foot->getRegion());

        double sumR = 0, sumG = 0, sumB = 0;  // sum of all non-saturated adjoining pixels
        double maxR = 0, maxG = 0, maxB = 0;  // maximum of non-saturated adjoining pixels

        for (auto span = bigFoot->getSpans()->begin(), send = bigFoot->getSpans()->end(); span != send;
             ++span) {
            int const y = span->getY() - y0;
            if (y < 0 || y >= height) {
                continue;
            }
            int sx0 = span->getX0() - x0;
            if (sx0 < 0) {
                sx0 = 0;
            }
            int sx1 = span->getX1() - x0;
            if (sx1 >= width) {
                sx1 = width - 1;
            }

            for (typename ImageT::iterator rptr = rim.at(sx0, y), rend = rim.at(sx1 + 1, y),
                                           gptr = gim.at(sx0, y), bptr = bim.at(sx0, y);
                 rptr != rend; ++rptr, ++gptr, ++bptr) {
                if (!((rptr.mask() | gptr.mask() | bptr.mask()) & SAT)) {
                    float val = rptr.image();
                    sumR += val;
                    if (val > maxR) {
                        maxR = val;
                    }

                    val = gptr.image();
                    sumG += val;
                    if (val > maxG) {
                        maxG = val;
                    }

                    val = bptr.image();
                    sumB += val;
                    if (val > maxB) {
                        maxB = val;
                    }
                }
            }
        }
        // OK, we have the mean fluxes for the pixels surrounding this set of saturated pixels
        // so we can figure out the proper values to use for the saturated ones
        float R = 0, G = 0, B = 0;  // mean intensities
        if (sumR + sumB + sumG > 0) {
            if (sumR > sumG) {
                if (sumR > sumB) {
                    R = useMaxPixel ? maxR : saturatedPixelValue;

                    G = (R * sumG) / sumR;
                    B = (R * sumB) / sumR;
                } else {
                    B = useMaxPixel ? maxB : saturatedPixelValue;
                    R = (B * sumR) / sumB;
                    G = (B * sumG) / sumB;
                }
            } else {
                if (sumG > sumB) {
                    G = useMaxPixel ? maxG : saturatedPixelValue;
                    R = (G * sumR) / sumG;
                    B = (G * sumB) / sumG;
                } else {
                    B = useMaxPixel ? maxB : saturatedPixelValue;
                    R = (B * sumR) / sumB;
                    G = (B * sumG) / sumB;
                }
            }
        }
        // Now that we know R, G, and B we can fix the values
        setR.setValue(R);
        foot->getSpans()->applyFunctor(setR, *rim.getImage());
        setG.setValue(G);
        foot->getSpans()->applyFunctor(setG, *gim.getImage());
        setB.setValue(B);
        foot->getSpans()->applyFunctor(setB, *bim.getImage());
    }
}

template void replaceSaturatedPixels(image::MaskedImage<float>& rim, image::MaskedImage<float>& gim,
                                     image::MaskedImage<float>& bim, int borderWidth,
                                     float saturatedPixelValue);
}
}
}