35 #include "boost/format.hpp"
42 namespace ex = pex::exceptions;
47 : _style(style), _orderX(
orderX), _orderY(orderY < 0 ?
orderX : orderY), _weighting(weighting) {
48 if (_orderX != _orderY) {
51 "due to a limitation in math::Chebyshev1Function2") %
60 template <
typename PixelT>
61 class ApproximateChebyshev :
public Approximate<PixelT> {
70 ~ApproximateChebyshev()
override;
79 int orderX,
int orderY)
const override;
81 int orderX,
int orderY)
const override;
87 inline void solveMatrix_Eigen(Eigen::MatrixXd&
a, Eigen::VectorXd&
b, Eigen::Map<Eigen::VectorXd>& c) {
88 Eigen::PartialPivLU<Eigen::MatrixXd> lu(
a);
96 template <
typename PixelT>
97 ApproximateChebyshev<PixelT>::ApproximateChebyshev(
102 ApproximateControl
const& ctrl
104 : Approximate<
PixelT>(xVec, yVec,
bbox, ctrl),
105 _poly(math::Chebyshev1Function2<double>(ctrl.getOrderX(),
lsst::
geom::Box2D(
bbox))) {
113 int const nTerm = _poly.getNParameters();
127 for (
int i = 0; i != nTerm; ++i) {
128 termCoeffs[i].reserve(nData);
131 for (
int iy = 0; iy != im.
getHeight(); ++iy) {
132 double const y = yVec[iy];
134 for (
int ix = 0; ix != im.
getWidth(); ++ix) {
135 double const x = xVec[ix];
137 for (
int i = 0; i != nTerm; ++i) {
138 _poly.setParameter(i, 1.0);
139 termCoeffs[i].push_back(_poly(
x,
y));
140 _poly.setParameter(i, 0.0);
146 A.setZero(nTerm, nTerm);
153 for (
int iy = 0; iy != im.
getHeight(); ++iy) {
157 double const val =
ptr.image();
158 double const ivar = ctrl.getWeighting() ? 1 /
ptr.variance() : 1.0;
163 for (
int i = 0; i != nTerm; ++i) {
164 double const c_i = termCoeffs[i][alpha];
165 double const tmp = c_i * ivar;
167 A(i, i) += c_i * tmp;
168 for (
int j = 0; j < i; ++j) {
169 double const c_j = termCoeffs[j][alpha];
170 A(i, j) += c_j * tmp;
176 if (ctrl.getWeighting()) {
178 "No valid points to fit. Variance is likely zero. Try weighting=False");
181 "No valid points to fit (even though weighting is False). "
182 "Check that binSize & approxOrderX settings are appropriate for image size.");
186 for (
int j = 0; j != nTerm; ++j) {
187 for (
int i = j + 1; i != nTerm; ++i) {
195 Eigen::Map<Eigen::VectorXd> c(&cvec[0], nTerm);
197 solveMatrix_Eigen(A,
b, c);
199 _poly.setParameters(cvec);
203 template <
typename PixelT>
204 ApproximateChebyshev<PixelT>::~ApproximateChebyshev() =
default;
216 template <
typename PixelT>
218 ApproximateChebyshev<PixelT>::doGetImage(
int orderX,
int orderY)
const {
230 for (
int iy = 0; iy != im->
getHeight(); ++iy) {
254 template <
typename PixelT>
256 ApproximateChebyshev<PixelT>::doGetMaskedImage(
int orderX,
int orderY)
const {
262 for (
int iy = 0; iy != im->getHeight(); ++iy) {
266 for (
typename MImageT::Image::x_iterator
ptr = im->row_begin(iy),
end = im->row_end(iy);
ptr !=
end;
278 template <
typename PixelT>
287 new ApproximateChebyshev<PixelT>(
x,
y, im,
bbox, ctrl));
298 #define INSTANTIATE(PIXEL_T) \
299 template std::shared_ptr<Approximate<PIXEL_T>> makeApproximate( \
300 std::vector<double> const& x, std::vector<double> const& y, \
301 image::MaskedImage<PIXEL_T> const& im, lsst::geom::Box2I const& bbox, \
302 ApproximateControl const& ctrl)
#define INSTANTIATE(FROMSYS, TOSYS)
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
A class to represent a 2-dimensional array of pixels.
A class to manipulate images, masks, and variance as a single object.
int getHeight() const
Return the number of rows in the image.
int getWidth() const
Return the number of columns in the image.
x_iterator row_end(int y) const
Return an x_iterator to the end of the image.
x_iterator row_begin(int y) const
Return an x_iterator to the start of the image.
ImagePtr getImage() const
Return a (shared_ptr to) the MaskedImage's image.
Control how to make an approximation.
int getOrderX() const
Return the order of approximation to use in the x-direction.
ApproximateControl(Style style, int orderX, int orderY=-1, bool weighting=true)
ctor
Style getStyle() const
Return the Style.
Style
Choose the type of approximation to use.
@ CHEBYSHEV
Use a 2-D Chebyshev polynomial.
int getOrderY() const
Return the order of approximation to use in the y-direction.
Approximate values for a MaskedImage.
ApproximateControl const _ctrl
desired approximation algorithm
lsst::geom::Box2I const _bbox
Domain for approximation.
2-dimensional weighted sum of Chebyshev polynomials of the first kind.
An integer coordinate rectangle.
Reports invalid arguments.
std::shared_ptr< Approximate< PixelT > > makeApproximate(std::vector< double > const &x, std::vector< double > const &y, image::MaskedImage< PixelT > const &im, lsst::geom::Box2I const &bbox, ApproximateControl const &ctrl)
Construct a new Approximate object, inferring the type from the type of the given MaskedImage.
Low-level polynomials (including special polynomials) in C++.
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
A base class for image defects.