PixelFlags.cc

Go to the documentation of this file.

// -*- lsst-c++ -*-
/*
 * LSST Data Management System
 * Copyright 2008-2013 LSST Corporation.
 *
 * This product includes software developed by the
 * LSST Project (http://www.lsst.org/).
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the LSST License Statement and
 * the GNU General Public License along with this program.  If not,
 * see <http://www.lsstcorp.org/LegalNotices/>.
 */

#include <cctype> // ::tolower
#include <algorithm> // std::transform
#include <cmath>

#include "ndarray/eigen.h"

#include "lsst/afw/detection/Psf.h"
#include "lsst/afw/detection/FootprintFunctor.h"
#include "lsst/afw/table/Source.h"
#include "lsst/meas/base/PixelFlags.h"

namespace lsst { namespace meas { namespace base {
namespace {
template <typename MaskedImageT>
class FootprintBits : public afw::detection::FootprintFunctor<MaskedImageT> {
public:
    explicit FootprintBits(MaskedImageT const& mimage) :
        afw::detection::FootprintFunctor<MaskedImageT>(mimage), _bits(0)
    {}

    void reset() {
        _bits = 0x0;
    }
    virtual void reset(afw::detection::Footprint const&) {}

    void operator()(typename MaskedImageT::xy_locator loc, 
                    int x,                                 
                    int y                                  
                   ) {
        _bits |= loc.mask(0, 0);
    }

    typename MaskedImageT::Mask::Pixel getBits() const { return _bits; }
private:
    typename MaskedImageT::Mask::Pixel _bits;
};

typedef afw::image::MaskedImage<float> MaskedImageF;

void updateFlags(PixelFlagsAlgorithm::KeyMap const & maskFlagToPixelFlag, const FootprintBits<MaskedImageF> & func,
                 afw::table::SourceRecord & measRecord) {
        for (auto const & i: maskFlagToPixelFlag) {
            try {
                if (func.getBits() & MaskedImageF::Mask::getPlaneBitMask(i.first)) {
                    measRecord.set(i.second, true);
                }
            }
            catch (pex::exceptions::InvalidParameterError & err) {
                throw LSST_EXCEPT(FatalAlgorithmError, err.what());
            }
        }
    }
} // end anonymous namespace

PixelFlagsAlgorithm::PixelFlagsAlgorithm(
    Control const & ctrl,
    std::string const & name,
    afw::table::Schema & schema
) : _ctrl(ctrl) {
    // Add generic keys first, which don't correspond to specific mask planes
    _generalFailureKey = schema.addField<afw::table::Flag>(name + "_flag",
                                        "general failure flag, set if anything went wring");
    _offImageKey = schema.addField<afw::table::Flag>(name+"_flag" + "_offimage",
                                        "Source center is off image");
    // Set all the flags that correspond to mask planes anywhere in the footprint
    _anyKeys["EDGE"] = schema.addField<afw::table::Flag>(name + "_flag_edge",
                                        "Source is outside usable exposure region (masked EDGE or NO_DATA)");
    _anyKeys["INTRP"] = schema.addField<afw::table::Flag>(name + "_flag_interpolated",
                                        "Interpolated pixel in the Source footprint");
    _anyKeys["SAT"] = schema.addField<afw::table::Flag>(name + "_flag_saturated",
                                        "Saturated pixel in the Source footprint");
    _anyKeys["CR"] = schema.addField<afw::table::Flag>(name + "_flag_cr",
                                        "Cosmic ray in the Source footprint");
    _anyKeys["BAD"] = schema.addField<afw::table::Flag>(name + "_flag_bad",
                                        "Bad pixel in the Source footprint");
    _anyKeys["SUSPECT"] = schema.addField<afw::table::Flag>(name + "_flag_suspect",
                                        "Source's footprint includes suspect pixels");
    // Flags that correspond to mask bits which occur in the center of the object
    _centerKeys["INTRP"] = schema.addField<afw::table::Flag>(name + "_flag_interpolatedCenter",
                                        "Interpolated pixel in the Source center");
    _centerKeys["SAT"] = schema.addField<afw::table::Flag>(name + "_flag_saturatedCenter",
                                        "Saturated pixel in the Source center");
    _centerKeys["CR"] = schema.addField<afw::table::Flag>(name + "_flag_crCenter",
                                        "Cosmic ray in the Source center");
    _centerKeys["SUSPECT"] = schema.addField<afw::table::Flag>(name + "_flag_suspectCenter",
                                        "Source's center is close to suspect pixels");

    // Read in the flags passed from the configuration, and add them to the schema
    for (auto const & i: _ctrl.masksFpCenter) {
        std::string maskName(i);
        std::transform(maskName.begin(), maskName.end(), maskName.begin(), ::tolower);
        _centerKeys[i] = schema.addField<afw::table::Flag>(name + "_flag_" + maskName + "Center",
                                        "Source center is close to "+ i + " pixels");
    }
    
    for (auto const & i: _ctrl.masksFpAnywhere) {
        std::string maskName(i);
        std::transform(maskName.begin(), maskName.end(), maskName.begin(), ::tolower);
        _anyKeys[i] = schema.addField<afw::table::Flag>(name + "_flag_" + maskName,
                                        "Source footprint includes " + i + " pixels");
    }
}

void PixelFlagsAlgorithm::measure(
    afw::table::SourceRecord & measRecord,
    afw::image::Exposure<float> const & exposure
) const {

    MaskedImageF mimage = exposure.getMaskedImage();
    FootprintBits<MaskedImageF> func(mimage);

    // Check if the measRecord has a valid centroid key, i.e. it was centroided
    afw::geom::Point2D center;
    if (measRecord.getTable()->getCentroidKey().isValid()) {
        center = measRecord.getCentroid();
        //  Catch NAN in centroid estimate
        if (std::isnan(center.getX()) || std::isnan(center.getY())) {
            throw LSST_EXCEPT(pex::exceptions::RuntimeError,
                              "Center point passed to PixelFlagsAlgorithm is NaN");
        }
    }
    else {
        // Set the general failure flag because using the Peak might affect
        // the current measurement
        measRecord.set(_generalFailureKey, true);
        // Attempting to set pixel flags with no centroider, the center will
        // be determined though the peak pixel. The first peak in the
        // footprint (supplied by the measurement framework) should be the
        // highest peak so that one will be used as a proxy for the central
        // tendency of the distribution of flux for the record.
        PTR(afw::detection::Footprint) footprint = measRecord.getFootprint();
        // If there is no footprint or the footprint contains no peaks, throw
        // a runtime error.
        if (!footprint || footprint->getPeaks().empty()) {
            throw LSST_EXCEPT(pex::exceptions::RuntimeError, "No footprint, or no footprint peaks detected");
        }
        else {
            center.setX(footprint->getPeaks().front().getFx());
            center.setY(footprint->getPeaks().front().getFy());
        }
    }

    //  Catch centroids off the image
    if (!mimage.getBBox().contains(afw::geom::Point2I(center))) {
        measRecord.set(_offImageKey, true);
        measRecord.set(_anyKeys.at("EDGE"), true);
    }

    // Check for bits set in the source's Footprint
    afw::detection::Footprint const & footprint(*measRecord.getFootprint());
    func.apply(footprint);

    // Set the EDGE flag if the bitmask has NO_DATA set
    try {
        if (func.getBits() & MaskedImageF::Mask::getPlaneBitMask("NO_DATA")) {
            measRecord.set(_anyKeys.at("EDGE"), true);
        }
    } catch (pex::exceptions::InvalidParameterError & err) {
        throw LSST_EXCEPT(FatalAlgorithmError, err.what());
    }

    // update the source record for the any keys
    updateFlags(_anyKeys, func, measRecord);

    // Check for bits set in the 3x3 box around the center
    afw::geom::Point2I llc(afw::image::positionToIndex(center.getX()) - 1,
                           afw::image::positionToIndex(center.getY()) - 1);

    afw::detection::Footprint const middle(afw::geom::BoxI(llc, afw::geom::ExtentI(3))); // central 3x3
    func.apply(middle);

    // Update the flags which have to do with the center of the footprint
    updateFlags(_centerKeys, func, measRecord);
}

void PixelFlagsAlgorithm::fail(afw::table::SourceRecord & measRecord, MeasurementError * error) const {
    measRecord.set(_generalFailureKey, true);
}

}}} // namespace lsst::meas::base