ellipse.h

// -*- LSST-C++ -*-
/*
 * This file is part of gauss2d.
 *
 * Developed for the LSST Data Management System.
 * This product includes software developed by the LSST Project
 * (https://www.lsst.org).
 * See the COPYRIGHT file at the top-level directory of this distribution
 * for details of code ownership.
 *
 * 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 GNU General Public License
 * along with this program.  If not, see <https://www.gnu.org/licenses/>.
 */
 
#ifndef LSST_GAUSS2D_ELLIPSE_H
#define LSST_GAUSS2D_ELLIPSE_H
 
#include <array>
#include <cmath>
#include <memory>
#include <stdexcept>
 
#include "object.h"
#include "to_string.h"
 
namespace lsst::gauss2d {
 
const double M_HWHM_SIGMA = 1.1774100225154746910115693264599;
const double M_SIGMA_HWHM = 0.84932180028801904272150283410295;
const double M_PI_180 = M_PI / 180.;
const double M_180_PI = 180. / M_PI;
 
class Ellipse;
class EllipseMajor;

class Covariance : public Object {
public:
    explicit Covariance(double sigma_x_sq = 0, double sigma_y_sq = 0, double cov_xy = 0);
    explicit Covariance(const Ellipse& ell);

    static void check(double sigma_x_sq, double sigma_y_sq, double cov_xy);
    void convolve(const Covariance& cov);
    double get_sigma_x_sq() const { return _sigma_x_sq; };
    double get_sigma_y_sq() const { return _sigma_y_sq; };
    double get_cov_xy() const { return _cov_xy; };
    std::array<double, 3> get_xyc() const { return {_sigma_x_sq, _sigma_y_sq, _cov_xy}; }

    std::shared_ptr<Covariance> make_convolution(const Covariance& cov) const;
    std::unique_ptr<Covariance> make_convolution_uniq(const Covariance& cov) const;

    void set(const Ellipse& ellipse);
    void set(double sigma_x_sq = 0, double sigma_y_sq = 0, double cov_xy = 0);
    void set_sigma_x_sq(double sigma_x_sq);
    void set_sigma_y_sq(double sigma_y_sq);
    void set_cov_xy(double cov_xy);
    void set_xyc(const std::array<double, 3>& xyc);
 
    std::string repr(bool name_keywords = false,
                     std::string_view namespace_separator = Object::CC_NAMESPACE_SEPARATOR) const override;
    std::string str() const override;
 
    bool operator==(const Covariance& other) const;
    bool operator!=(const Covariance& other) const;
 
    friend std::ostream& operator<<(std::ostream& out, const Covariance& obj);
 
private:
    double _sigma_x_sq = 0;
    double _sigma_y_sq = 0;
    double _cov_xy = 0;
};

class EllipseData : public Object {
public:
    virtual ~EllipseData() = default;

    static void check(double size_x, double size_y, double rho, std::string_view error_suffix = "") {
        if (!(size_x >= 0) || !(size_y >= 0) || !(rho >= -1 && rho <= 1)) {
            throw std::invalid_argument("Invalid size_x, size_y, rho=" + to_string_float(size_x) + ","
                                        + to_string_float(size_y) + "," + to_string_float(rho)
                                        + "; sigma_x,y >= 0 and 1 >= rho >= -1 required."
                                        + std::string(error_suffix));
        }
    }

    static void check_size(double size, std::string_view error_suffix = "") {
        if (!(size >= 0)) {
            throw std::invalid_argument("Invalid size=" + to_string_float(size) + "; size >= 0 required."
                                        + std::string(error_suffix));
        }
    }

    static void check_rho(double rho, std::string_view error_suffix = "") {
        if (!(rho >= -1 && rho <= 1)) {
            throw std::invalid_argument("Invalid rho=" + to_string_float(rho) + "; 1 >= rho >= -1 required."
                                        + std::string(error_suffix));
        }
    }

    virtual void convolve(const Ellipse& ell);
    virtual double get_area() const;
    virtual double get_cov_xy() const;
    virtual double get_hwhm_x() const;
    virtual double get_hwhm_y() const;
    virtual double get_radius_trace() const;
    virtual double get_sigma_x() const = 0;
    virtual double get_sigma_y() const = 0;
    virtual double get_rho() const = 0;
    virtual double get_sigma_x_sq() const;
    virtual double get_sigma_y_sq() const;
    virtual double get_sigma_xy() const;
    virtual std::array<double, 3> get_hxyr() const;
    virtual std::array<double, 3> get_xyr() const;

    virtual void set(double sigma_x, double sigma_y, double rho);
    virtual void set(const Covariance& covar);
    virtual void set(const EllipseMajor& ellipse);
    virtual void set_h(double hwhm_x, double hwhm_y, double rho);
    virtual void set_hwhm_x(double hwhm_x);
    virtual void set_hwhm_y(double hwhm_y);
    virtual void set_sigma_x(double sigma_x) = 0;
    virtual void set_sigma_y(double sigma_y) = 0;
    virtual void set_rho(double rho) = 0;
    virtual void set_hxyr(const std::array<double, 3>& hxyr);
    virtual void set_xyr(const std::array<double, 3>& xyr);
 
    virtual std::string repr(bool name_keywords = false, std::string_view namespace_separator
                                                         = Object::CC_NAMESPACE_SEPARATOR) const override
            = 0;
    virtual std::string str() const override = 0;
 
    bool operator==(const EllipseData& other) const { return get_xyr() == other.get_xyr(); };
    bool operator!=(const EllipseData& other) const { return !(*this == other); };
 
    friend std::ostream& operator<<(std::ostream& out, const EllipseData& obj) {
        out << obj.str();
        return out;
    }
};

class EllipseValues : public EllipseData {
public:
    explicit EllipseValues(std::shared_ptr<double> sigma_x, std::shared_ptr<double> sigma_y,
                           std::shared_ptr<double> rho = nullptr)
            : _sigma_x(sigma_x == nullptr ? std::make_shared<double>(0) : std::move(sigma_x)),
              _sigma_y(sigma_y == nullptr ? std::make_shared<double>(0) : std::move(sigma_y)),
              _rho(rho == nullptr ? std::make_shared<double>(0) : std::move(rho)) {};

    explicit EllipseValues(double sigma_x = 0, double sigma_y = 0, double rho = 0)
            : _sigma_x(std::make_shared<double>(sigma_x)),
              _sigma_y(std::make_shared<double>(sigma_y)),
              _rho(std::make_shared<double>(rho)) {};
 
    double get_sigma_x() const override;
    double get_sigma_y() const override;
    double get_rho() const override;
    std::array<double, 3> get_xyr() const override;
 
    void set(double sigma_x, double sigma_y, double rho) override;
    void set_h(double hwhm_x, double hwhm_y, double rho) override;
    void set_sigma_x(double sigma_x) override;
    void set_sigma_y(double sigma_y) override;
    void set_rho(double rho) override;
 
    std::string repr(bool name_keywords = false,
                     std::string_view namespace_separator = Object::CC_NAMESPACE_SEPARATOR) const override;
    std::string str() const override;
 
private:
    std::shared_ptr<double> _sigma_x;
    std::shared_ptr<double> _sigma_y;
    std::shared_ptr<double> _rho;
};

class Ellipse : public EllipseData {
public:
    explicit Ellipse(std::shared_ptr<EllipseData> data);
    explicit Ellipse(double sigma_x = 0, double sigma_y = 0, double rho = 0);
    explicit Ellipse(const Covariance& covar);
    explicit Ellipse(const EllipseMajor& ellipse);

    const EllipseData& get_data() const;
    double get_rho() const override;
    double get_sigma_x() const override;
    double get_sigma_y() const override;
    std::shared_ptr<Ellipse> make_convolution(const Ellipse& ell) const;
    std::unique_ptr<Ellipse> make_convolution_uniq(const Ellipse& ell) const;
 
    void set_rho(double rho) override;
    void set_sigma_x(double sigma_x) override;
    void set_sigma_y(double sigma_y) override;
 
    std::string repr(bool name_keywords = false,
                     std::string_view namespace_separator = Object::CC_NAMESPACE_SEPARATOR) const override;
    std::string str() const override;
 
    bool operator==(const Ellipse& other) const;
    bool operator!=(const Ellipse& other) const;
 
private:
    std::shared_ptr<EllipseData> _data;
};

class EllipseMajor : public Object {
public:
    explicit EllipseMajor(double r_major, double axrat, double angle, bool degrees = false);
    explicit EllipseMajor(const Covariance& covar, bool degrees = false);
    explicit EllipseMajor(const Ellipse& ellipse, bool degrees = false);

    static void check(double r_major, double axrat, double angle) {
        if (!(r_major >= 0) || !(axrat >= 0 && axrat <= 1)) {
            throw std::invalid_argument("Invalid r_major, axrat, angle=" + std::to_string(r_major) + ","
                                        + std::to_string(axrat) + "," + std::to_string(angle)
                                        + "; r_major >= 0, 1 >= axrat >= 0 required.");
        }
    }

    double get_area() const { return M_PI * _r_major * _r_major * _axrat; }
    double get_r_major() const { return _r_major; }
    double get_axrat() const { return _axrat; }
    double get_angle() const { return _angle; }
    double get_angle_degrees() const { return _degrees ? _angle : _angle * M_180_PI; }
    double get_angle_radians() const { return _degrees ? _angle * M_PI_180 : _angle; }
    std::array<double, 3> get_rqa() const { return {_r_major, _axrat, _angle}; }
    bool is_degrees() const { return _degrees; }
 
    void set(double r_major, double axrat, double angle);
    void set_r_major(double r_major);
    void set_axrat(double axrat);
    void set_angle(double angle);
    void set_degrees(bool degrees);
    void set_rqa(const std::array<double, 3>& rqa);
 
    std::string repr(bool name_keywords = false,
                     std::string_view namespace_separator = Object::CC_NAMESPACE_SEPARATOR) const override;
    std::string str() const override;
 
    bool operator==(const EllipseMajor& other) const;
    bool operator!=(const EllipseMajor& other) const;
 
private:
    double _r_major = 0.;
    double _axrat = 1.;
    double _angle = 0.;
    bool _degrees = false;
};
 
}  // namespace lsst::gauss2d
#endif