SdssShape.h

Go to the documentation of this file.

// -*- lsst-c++ -*-
/*
 * LSST Data Management System
 * Copyright 2008-2016 AURA/LSST.
 *
 * 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 <https://www.lsstcorp.org/LegalNotices/>.
 */

#ifndef LSST_MEAS_BASE_SdssShape_h_INCLUDED
#define LSST_MEAS_BASE_SdssShape_h_INCLUDED

#include <bitset>

#include "lsst/pex/config.h"
#include "lsst/afw/geom/SkyWcs.h"
#include "lsst/afw/image/PhotoCalib.h"
#include "lsst/afw/image/Exposure.h"
#include "lsst/afw/geom/ellipses/Quadrupole.h"
#include "lsst/afw/table/aggregates.h"
#include "lsst/meas/base/FluxUtilities.h"
#include "lsst/meas/base/CentroidUtilities.h"
#include "lsst/meas/base/ShapeUtilities.h"
#include "lsst/meas/base/InputUtilities.h"
#include "lsst/meas/base/Algorithm.h"

namespace lsst {
namespace meas {
namespace base {

class SdssShapeResult;

class SdssShapeControl {
public:
    LSST_CONTROL_FIELD(background, double, "Additional value to add to background");
    LSST_CONTROL_FIELD(maxIter, int, "Maximum number of iterations");
    LSST_CONTROL_FIELD(maxShift, double, "Maximum centroid shift, limited to 2-10");
    LSST_CONTROL_FIELD(tol1, float, "Convergence tolerance for e1,e2");
    LSST_CONTROL_FIELD(tol2, float, "Convergence tolerance for FWHM");
    LSST_CONTROL_FIELD(doMeasurePsf, bool, "Whether to also compute the shape of the PSF model");

    SdssShapeControl()
            : background(0.0), maxIter(100), maxShift(), tol1(1E-5), tol2(1E-4), doMeasurePsf(true) {}
};

class SdssShapeResultKey : public afw::table::FunctorKey<SdssShapeResult> {
public:
    static SdssShapeResultKey addFields(afw::table::Schema& schema, std::string const& name,
                                        bool doMeasurePsf);

    SdssShapeResultKey() {}

    SdssShapeResultKey(afw::table::SubSchema const& s);

    virtual SdssShapeResult get(afw::table::BaseRecord const& record) const;

    virtual void set(afw::table::BaseRecord& record, SdssShapeResult const& value) const;

    virtual afw::geom::ellipses::Quadrupole getPsfShape(afw::table::BaseRecord const& record) const;

    virtual void setPsfShape(afw::table::BaseRecord& record,
                             afw::geom::ellipses::Quadrupole const& value) const;

    bool operator==(SdssShapeResultKey const& other) const;
    bool operator!=(SdssShapeResultKey const& other) const { return !(*this == other); }

    bool isValid() const;

    FlagHandler const& getFlagHandler() const { return _flagHandler; }

private:
    bool _includePsf;
    ShapeResultKey _shapeResult;
    CentroidResultKey _centroidResult;
    FluxResultKey _instFluxResult;
    afw::table::QuadrupoleKey _psfShapeResult;
    afw::table::Key<ErrElement> _instFlux_xx_Cov;
    afw::table::Key<ErrElement> _instFlux_yy_Cov;
    afw::table::Key<ErrElement> _instFlux_xy_Cov;
    FlagHandler _flagHandler;
};

class SdssShapeAlgorithm : public SimpleAlgorithm {
public:
    // Structures and routines to manage flaghandler
    static FlagDefinitionList const& getFlagDefinitions();
    static unsigned int const N_FLAGS = 6;
    static FlagDefinition const FAILURE;
    static FlagDefinition const UNWEIGHTED_BAD;
    static FlagDefinition const UNWEIGHTED;
    static FlagDefinition const SHIFT;
    static FlagDefinition const MAXITER;
    static FlagDefinition const PSF_SHAPE_BAD;

    typedef SdssShapeControl Control;
    typedef SdssShapeResult Result;
    typedef SdssShapeResultKey ResultKey;

    // NOTE: In order to accommodate the optional setting of additional fields when running with
    //       doMeasurePsf = true (do set extra fields) or false (do NOT set extra fields), all of
    //       the code in SdssShape assumes that PSF_SHAPE_BAD is the last entry in the enum list.
    //       If new flags are added, be sure to add them above the PSF_SHAPE_BAD entry.

    SdssShapeAlgorithm(Control const& ctrl, std::string const& name, afw::table::Schema& schema);

    template <typename ImageT>
    static Result computeAdaptiveMoments(ImageT const& image, geom::Point2D const& position,
                                         bool negative = false, Control const& ctrl = Control());

    template <typename ImageT>
    static FluxResult computeFixedMomentsFlux(ImageT const& image,
                                              afw::geom::ellipses::Quadrupole const& shape,
                                              geom::Point2D const& position);

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

    virtual void fail(afw::table::SourceRecord& measRecord, MeasurementError* error = nullptr) const;

private:
    Control _ctrl;
    ResultKey _resultKey;
    SafeCentroidExtractor _centroidExtractor;
};

class SdssShapeResult : public ShapeResult, public CentroidResult, public FluxResult {
public:
    ErrElement instFlux_xx_Cov;  
    ErrElement instFlux_yy_Cov;  
    ErrElement instFlux_xy_Cov;  

    std::bitset<SdssShapeAlgorithm::N_FLAGS> flags;  

    // TODO is this workaround still needed?
    bool getFlag(unsigned int index) const { return flags[index]; }

    bool getFlag(std::string const& name) const {
        return flags[SdssShapeAlgorithm::getFlagDefinitions().getDefinition(name).number];
    }

    SdssShapeResult();  
};

class SdssShapeTransform : public BaseTransform {
public:
    typedef SdssShapeControl Control;

    SdssShapeTransform(Control const& ctrl, std::string const& name, afw::table::SchemaMapper& mapper);

    /*
     * @brief Perform transformation from inputCatalog to outputCatalog.
     *
     * @param[in]     inputCatalog   Source of data to be transformed
     * @param[in,out] outputCatalog  Container for transformed results
     * @param[in]     wcs            World coordinate system under which transformation will take place
     * @param[in]     photoCalib     Photometric calibration under which transformation will take place
     * @throws        LengthError    Catalog sizes do not match
     */
    virtual void operator()(afw::table::SourceCatalog const& inputCatalog,
                            afw::table::BaseCatalog& outputCatalog, afw::geom::SkyWcs const& wcs,
                            afw::image::PhotoCalib const& photoCalib) const;

private:
    FluxTransform _instFluxTransform;
    CentroidTransform _centroidTransform;
    ShapeResultKey _outShapeKey;
    afw::table::QuadrupoleKey _outPsfShapeKey;
    bool _transformPsf;
};

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

#endif  // !LSST_MEAS_BASE_SdssShape_h_INCLUDED