Skip to content

Commit

Permalink
Change polynomial degree to `order.
Browse files Browse the repository at this point in the history
After a straw poll on slack and subsequent discussion, I've made jointcal
consistent with the rest of the stack by changing `degree` to `order` as the
term for the maximum monomial.
  • Loading branch information
parejkoj committed Apr 17, 2018
1 parent c82a162 commit 430e93e
Show file tree
Hide file tree
Showing 21 changed files with 237 additions and 239 deletions.
4 changes: 2 additions & 2 deletions include/lsst/jointcal/ConstrainedAstrometryModel.h
Original file line number Diff line number Diff line change
Expand Up @@ -31,8 +31,8 @@ namespace jointcal {
class ConstrainedAstrometryModel : public AstrometryModel {
public:
ConstrainedAstrometryModel(CcdImageList const &ccdImageList,
std::shared_ptr<ProjectionHandler const> projectionHandler, int chipDegree,
int visitDegree);
std::shared_ptr<ProjectionHandler const> projectionHandler, int chipOrder,
int visitOrder);

/// No copy or move: there is only ever one instance of a given model (i.e. per ccd+visit)
ConstrainedAstrometryModel(ConstrainedAstrometryModel const &) = delete;
Expand Down
4 changes: 2 additions & 2 deletions include/lsst/jointcal/ConstrainedPhotometryModel.h
Original file line number Diff line number Diff line change
Expand Up @@ -30,10 +30,10 @@ class ConstrainedPhotometryModel : public PhotometryModel {
* @param ccdImageList The list of CCDImages to construct the model for.
* @param focalPlaneBBox The bounding box of the camera's focal plane, defining the domain of the
* visit polynomial.
* @param[in] visitDegree The degree of the visit polynomial.
* @param[in] visitOrder The order of the visit polynomial.
*/
explicit ConstrainedPhotometryModel(CcdImageList const &ccdImageList,
afw::geom::Box2D const &focalPlaneBBox, int visitDegree = 7);
afw::geom::Box2D const &focalPlaneBBox, int visitOrder = 7);

/// No copy or move: there is only ever one instance of a given model (i.e. per ccd+visit)
ConstrainedPhotometryModel(ConstrainedPhotometryModel const &) = delete;
Expand Down
36 changes: 19 additions & 17 deletions include/lsst/jointcal/Gtransfo.h
Original file line number Diff line number Diff line change
Expand Up @@ -252,26 +252,31 @@ bool isIntegerShift(const Gtransfo *gtransfo);
//! Polynomial transformation class.
class GtransfoPoly : public Gtransfo {
public:
//! Default transfo : identity for all degrees (>=1 ). The degree refers to the highest total power (x+y)
//! of monomials.
GtransfoPoly(const unsigned degree = 1);
/**
* Default transfo : identity for all orders (>=1 ).
*
* @param order The highest total power (x+y) of monomials of this polynomial.
*/
GtransfoPoly(const unsigned order = 1);

//! Constructs a "polynomial image" from an existing transfo, over a specified domain
GtransfoPoly(const Gtransfo *gtransfo, const Frame &frame, unsigned degree, unsigned nPoint = 1000);
GtransfoPoly(const Gtransfo *gtransfo, const Frame &frame, unsigned order, unsigned nPoint = 1000);

/**
* Constructs a polynomial approximation to an afw::geom::TransformPoint2ToPoint2.
*
* @param[in] transform The transform to be approximated.
* @param[in] domain The valid domain of the transform.
* @param[in] degree The polynomial degree to use when approximating.
* @param[in] order The polynomial order to use when approximating.
* @param[in] nSteps The number of sample points per axis (nSteps^2 total points).
*/
GtransfoPoly(std::shared_ptr<afw::geom::TransformPoint2ToPoint2> transform, jointcal::Frame const &domain,
unsigned const degree, unsigned const nSteps = 50);
unsigned const order, unsigned const nSteps = 50);

// sets the polynomial degree.
void setDegree(const unsigned degree);
/// Sets the polynomial order (the highest sum of exponents of the largest monomial).
void setOrder(const unsigned order);
/// Returns the polynomial order.
unsigned getOrder() const { return _order; }

using Gtransfo::apply; // to unhide Gtransfo::apply(Point const &)

Expand All @@ -284,9 +289,6 @@ class GtransfoPoly : public Gtransfo {
//! a mix of apply and Derivative
virtual void transformPosAndErrors(const FatPoint &in, FatPoint &out) const override;

//! returns degree
unsigned getDegree() const { return _degree; }

//! total number of parameters
int getNpar() const override { return 2 * _nterms; }

Expand Down Expand Up @@ -319,7 +321,7 @@ class GtransfoPoly : public Gtransfo {
//! write access
double &coeff(const unsigned powX, const unsigned powY, const unsigned whichCoord);

//! read access, zero if beyond degree
//! read access, zero if beyond order
double coeffOrZero(const unsigned powX, const unsigned powY, const unsigned whichCoord) const;

double determinant() const;
Expand All @@ -343,7 +345,7 @@ class GtransfoPoly : public Gtransfo {
private:
double computeFit(StarMatchList const &starMatchList, Gtransfo const &InTransfo, const bool UseErrors);

unsigned _degree; // the degree
unsigned _order; // The highest sum of exponents of the largest monomial.
unsigned _nterms; // number of parameters per coordinate
std::vector<double> _coeffs; // the actual coefficients
// both polynomials in a single vector to speed up allocation and copies
Expand Down Expand Up @@ -373,13 +375,13 @@ class GtransfoPoly : public Gtransfo {
* @param forward Transform to be inverted.
* @param[in] domain The domain of forward.
* @param[in] precision Require that \f$chi2 / (nsteps^2) < precision^2\f$.
* @param[in] maxDegree The maximum degree allowed of the inverse polynomial.
* @param[in] maxOrder The maximum order allowed of the inverse polynomial.
* @param[in] nSteps The number of sample points per axis (nSteps^2 total points).
*
* @return A polynomial that best approximates forward.
*/
std::shared_ptr<GtransfoPoly> inversePolyTransfo(Gtransfo const &forward, Frame const &domain,
double const precision, int const maxDegree = 9,
double const precision, int const maxOrder = 9,
unsigned const nSteps = 50);

GtransfoLin normalizeCoordinatesTransfo(const Frame &frame);
Expand All @@ -393,7 +395,7 @@ class GtransfoLin : public GtransfoPoly {
//! the default constructor constructs the do-nothing transformation.
GtransfoLin() : GtransfoPoly(1){};

//! This triggers an exception if P.degree() != 1
//! This triggers an exception if P.getOrder() != 1
explicit GtransfoLin(GtransfoPoly const &gtransfoPoly);

//! enables to combine linear tranformations: T1=T2*T3 is legal.
Expand Down Expand Up @@ -444,7 +446,7 @@ class GtransfoLin : public GtransfoPoly {
friend class GtransfoPoly; // // for Gtransfo::Derivative

private:
void setDegree(const unsigned degree); // to hide GtransfoPoly::setDegree
void setOrder(const unsigned order); // to hide GtransfoPoly::setOrder
};

/*=============================================================*/
Expand Down
20 changes: 10 additions & 10 deletions include/lsst/jointcal/PhotometryTransfo.h
Original file line number Diff line number Diff line change
Expand Up @@ -129,36 +129,36 @@ class PhotometryTransfoSpatiallyInvariant : public PhotometryTransfo {
};

/**
* nth-degree 2d Chebyshev photometry transfo.
* nth-order 2d Chebyshev photometry transfo.
*
* The 2-d Chebyshev polynomial used here is defined as:
*
* @f[
* f(x,y) = \sum_i \sum_j a_{i,j} T_i(x) T_j(y)
* @f]
*
* where @f$T_n(x)@f$ is the n-th degree Chebyshev polynomial of @f$x@f$ and
* where @f$T_n(x)@f$ is the n-th order Chebyshev polynomial of @f$x@f$ and
* @f$a_{i,j}@f$ is the corresponding coefficient of the (i,j) polynomial term.
*
* Note that the polynomial degree=n means that the highest terms will be of the form:
* Note that the polynomial order=n means that the highest terms will be of the form:
* @f[
* a_{0,n}*x^n*y^0, a_{n-1,1}*x^(n-1)*y^1, ..., a_{1,n-1}*x^1*y^(n-1), a_{n,0}*x^0*y^n
* @f]
*/
class PhotometryTransfoChebyshev : public PhotometryTransfo {
public:
/**
* Create an identity (a_0,0==1) Chebyshev transfo with terms up to degree in (x*y).
* Create an identity (a_0,0==1) Chebyshev transfo with terms up to order in (x*y).
*
* @param[in] degree The maximum degree in (x*y).
* @param[in] order The maximum order in (x*y).
* @param[in] bbox The bounding box it is valid within, to rescale it to [-1,1].
*/
PhotometryTransfoChebyshev(size_t degree, afw::geom::Box2D const &bbox);
PhotometryTransfoChebyshev(size_t order, afw::geom::Box2D const &bbox);

/**
* Create a Chebyshev transfo with the specified coefficients.
*
* The polynomial degree is determined from the number of coefficients.
* The polynomial order is determined from the number of coefficients.
*
* @param coefficients The polynomial coefficients.
* @param[in] bbox The bounding box it is valid within, to rescale it to [-1,1].
Expand Down Expand Up @@ -193,7 +193,7 @@ class PhotometryTransfoChebyshev : public PhotometryTransfo {
/// @copydoc PhotometryTransfo::getParameters
Eigen::VectorXd getParameters() const override;

ndarray::Size getDegree() const { return _degree; }
ndarray::Size getOrder() const { return _order; }

afw::geom::Box2D getBBox() const { return _bbox; }

Expand All @@ -204,8 +204,8 @@ class PhotometryTransfoChebyshev : public PhotometryTransfo {
afw::geom::Box2D _bbox; // the domain of this function
afw::geom::AffineTransform _toChebyshevRange; // maps points from the bbox to [-1,1]x[-1,1]

ndarray::Array<double, 2, 2> _coefficients; // shape=(degree+1, degree+1)
ndarray::Size _degree;
ndarray::Array<double, 2, 2> _coefficients; // shape=(order+1, order+1)
ndarray::Size _order;
ndarray::Size _nParameters;

// Compute the integral of this function over its bounding-box.
Expand Down
2 changes: 1 addition & 1 deletion include/lsst/jointcal/SimpleAstrometryModel.h
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ class SimpleAstrometryModel : public AstrometryModel {
//! Sky2TP is just a name, it can be anything
SimpleAstrometryModel(CcdImageList const &ccdImageList,
const std::shared_ptr<ProjectionHandler const> projectionHandler, bool initFromWCS,
unsigned nNotFit = 0, unsigned degree = 3);
unsigned nNotFit = 0, unsigned order = 3);

/// No copy or move: there is only ever one instance of a given model (i.e.. per ccd+visit)
SimpleAstrometryModel(SimpleAstrometryModel const &) = delete;
Expand Down
4 changes: 2 additions & 2 deletions python/lsst/jointcal/astrometryModels.cc
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,7 @@ void declareSimpleAstrometryModel(py::module &mod) {

cls.def(py::init<CcdImageList const &, const std::shared_ptr<ProjectionHandler const>, bool, unsigned,
unsigned>(),
"ccdImageList"_a, "projectionHandler"_a, "initFromWcs"_a, "nNotFit"_a = 0, "degree"_a = 3);
"ccdImageList"_a, "projectionHandler"_a, "initFromWcs"_a, "nNotFit"_a = 0, "order"_a = 3);

cls.def("getTransfo", &SimpleAstrometryModel::getTransfo, py::return_value_policy::reference_internal);
}
Expand All @@ -67,7 +67,7 @@ void declareConstrainedAstrometryModel(py::module &mod) {
mod, "ConstrainedAstrometryModel");

cls.def(py::init<CcdImageList const &, std::shared_ptr<ProjectionHandler const>, int, int>(),
"ccdImageList"_a, "projectionHandler"_a, "chipDegree"_a, "visitDegree"_a);
"ccdImageList"_a, "projectionHandler"_a, "chipOrder"_a, "visitOrder"_a);

cls.def("getChipTransfo", &ConstrainedAstrometryModel::getChipTransfo,
py::return_value_policy::reference_internal);
Expand Down
6 changes: 3 additions & 3 deletions python/lsst/jointcal/gtransfo.cc
Original file line number Diff line number Diff line change
Expand Up @@ -61,8 +61,8 @@ void declareGtransfoIdentity(py::module &mod) {
void declareGtransfoPoly(py::module &mod) {
py::class_<GtransfoPoly, std::shared_ptr<GtransfoPoly>, Gtransfo> cls(mod, "GtransfoPoly");

cls.def(py::init<const unsigned>(), "degree"_a);
cls.def("getDegree", &GtransfoPoly::getDegree);
cls.def(py::init<const unsigned>(), "order"_a);
cls.def("getOrder", &GtransfoPoly::getOrder);
cls.def("coeff", (double (GtransfoPoly::*)(unsigned const, unsigned const, unsigned const) const) &
GtransfoPoly::coeff);
cls.def("write", [](GtransfoPoly const &self) {
Expand Down Expand Up @@ -139,7 +139,7 @@ PYBIND11_PLUGIN(gtransfo) {

// utility functions
mod.def("inversePolyTransfo", &inversePolyTransfo, "forward"_a, "domain"_a, "precision"_a,
"maxDegree"_a = 9, "nSteps"_a = 50);
"maxOrder"_a = 9, "nSteps"_a = 50);

return mod.ptr();
}
Expand Down
24 changes: 12 additions & 12 deletions python/lsst/jointcal/jointcal.py
Original file line number Diff line number Diff line change
Expand Up @@ -147,18 +147,18 @@ class JointcalConfig(pexConfig.Config):
dtype=int,
default=2,
)
astrometrySimpleDegree = pexConfig.Field(
doc="Polynomial degree for fitting the simple astrometry model.",
astrometrySimpleOrder = pexConfig.Field(
doc="Polynomial order for fitting the simple astrometry model.",
dtype=int,
default=3,
)
astrometryChipDegree = pexConfig.Field(
doc="Degree of the per-chip transform for the constrained astrometry model.",
astrometryChipOrder = pexConfig.Field(
doc="Order of the per-chip transform for the constrained astrometry model.",
dtype=int,
default=1,
)
astrometryVisitDegree = pexConfig.Field(
doc="Degree of the per-visit transform for the constrained astrometry model.",
astrometryVisitOrder = pexConfig.Field(
doc="Order of the per-visit transform for the constrained astrometry model.",
dtype=int,
default=5,
)
Expand All @@ -181,8 +181,8 @@ class JointcalConfig(pexConfig.Config):
allowed={"simple": "One constant zeropoint per ccd and visit",
"constrained": "Constrained zeropoint per ccd, and one polynomial per visit"}
)
photometryVisitDegree = pexConfig.Field(
doc="Degree of the per-visit polynomial transform for the constrained photometry model.",
photometryVisitOrder = pexConfig.Field(
doc="Order of the per-visit polynomial transform for the constrained photometry model.",
dtype=int,
default=7,
)
Expand Down Expand Up @@ -544,7 +544,7 @@ def _fit_photometry(self, associations, dataName=None):
if self.config.photometryModel == "constrained":
model = lsst.jointcal.ConstrainedPhotometryModel(associations.getCcdImageList(),
self.focalPlaneBBox,
visitDegree=self.config.photometryVisitDegree)
visitOrder=self.config.photometryVisitOrder)
elif self.config.photometryModel == "simple":
model = lsst.jointcal.SimplePhotometryModel(associations.getCcdImageList())

Expand Down Expand Up @@ -622,14 +622,14 @@ def _fit_astrometry(self, associations, dataName=None):
if self.config.astrometryModel == "constrained":
model = lsst.jointcal.ConstrainedAstrometryModel(associations.getCcdImageList(),
sky_to_tan_projection,
chipDegree=self.config.astrometryChipDegree,
visitDegree=self.config.astrometryVisitDegree)
chipOrder=self.config.astrometryChipOrder,
visitOrder=self.config.astrometryVisitOrder)
elif self.config.astrometryModel == "simple":
model = lsst.jointcal.SimpleAstrometryModel(associations.getCcdImageList(),
sky_to_tan_projection,
self.config.useInputWcs,
nNotFit=0,
degree=self.config.astrometrySimpleDegree)
order=self.config.astrometrySimpleOrder)

fit = lsst.jointcal.AstrometryFit(associations, model, self.config.posError)
chi2 = fit.computeChi2()
Expand Down
2 changes: 1 addition & 1 deletion python/lsst/jointcal/photometryModels.cc
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,7 @@ void declareConstrainedPhotometryModel(py::module &mod) {
py::class_<ConstrainedPhotometryModel, std::shared_ptr<ConstrainedPhotometryModel>, PhotometryModel> cls(
mod, "ConstrainedPhotometryModel");
cls.def(py::init<CcdImageList const &, afw::geom::Box2D const &, int>(), "CcdImageList"_a, "bbox"_a,
"visitDegree"_a = 7);
"visitOrder"_a = 7);
}

PYBIND11_PLUGIN(photometryModels) {
Expand Down
4 changes: 2 additions & 2 deletions python/lsst/jointcal/photometryTransfo.cc
Original file line number Diff line number Diff line change
Expand Up @@ -69,12 +69,12 @@ void declarePhotometryTransfoChebyshev(py::module &mod) {
py::class_<PhotometryTransfoChebyshev, std::shared_ptr<PhotometryTransfoChebyshev>, PhotometryTransfo>
cls(mod, "PhotometryTransfoChebyshev");

cls.def(py::init<size_t, afw::geom::Box2D const &>(), "degree"_a, "bbox"_a);
cls.def(py::init<size_t, afw::geom::Box2D const &>(), "order"_a, "bbox"_a);
cls.def(py::init<ndarray::Array<double, 2, 2> const &, afw::geom::Box2D const &>(), "coefficients"_a,
"bbox"_a);

cls.def("getCoefficients", &PhotometryTransfoChebyshev::getCoefficients);
cls.def("getDegree", &PhotometryTransfoChebyshev::getDegree);
cls.def("getOrder", &PhotometryTransfoChebyshev::getOrder);
cls.def("getBBox", &PhotometryTransfoChebyshev::getBBox);
}

Expand Down
8 changes: 4 additions & 4 deletions src/ConstrainedAstrometryModel.cc
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ namespace jointcal {

ConstrainedAstrometryModel::ConstrainedAstrometryModel(
CcdImageList const &ccdImageList, std::shared_ptr<ProjectionHandler const> projectionHandler,
int chipDegree, int visitDegree)
int chipOrder, int visitOrder)
: _sky2TP(projectionHandler) {
// keep track of which chip we want to hold fixed (the one closest to the middle of the focal plane)
double minRadius2 = std::numeric_limits<double>::infinity();
Expand All @@ -55,7 +55,7 @@ ConstrainedAstrometryModel::ConstrainedAstrometryModel(
auto chip = im.getCcdId();
auto visitp = _visitMap.find(visit);
if (visitp == _visitMap.end()) {
_visitMap[visit] = std::make_shared<SimplePolyMapping>(GtransfoLin(), GtransfoPoly(visitDegree));
_visitMap[visit] = std::make_shared<SimplePolyMapping>(GtransfoLin(), GtransfoPoly(visitOrder));
}
auto chipp = _chipMap.find(chip);
if (chipp == _chipMap.end()) {
Expand All @@ -70,7 +70,7 @@ ConstrainedAstrometryModel::ConstrainedAstrometryModel(
im.getDetector()->getTransform(afw::cameraGeom::PIXELS, afw::cameraGeom::FOCAL_PLANE);
Frame const &frame = im.getImageFrame();
// construct the chip gtransfo by approximating the pixel->Focal afw::geom::Transform.
GtransfoPoly pol = GtransfoPoly(pixelsToFocal, frame, chipDegree);
GtransfoPoly pol = GtransfoPoly(pixelsToFocal, frame, chipOrder);
GtransfoLin shiftAndNormalize = normalizeCoordinatesTransfo(frame);
_chipMap[chip] =
std::make_shared<SimplePolyMapping>(shiftAndNormalize, pol * shiftAndNormalize.invert());
Expand All @@ -91,7 +91,7 @@ ConstrainedAstrometryModel::ConstrainedAstrometryModel(
if (_chipMap.find(chip) == _chipMap.end()) {
LOGLS_WARN(_log, "Chip " << chip << " is missing in the reference exposure, expect troubles.");
GtransfoLin norm = normalizeCoordinatesTransfo(im.getImageFrame());
_chipMap[chip] = std::make_shared<SimplePolyMapping>(norm, GtransfoPoly(chipDegree));
_chipMap[chip] = std::make_shared<SimplePolyMapping>(norm, GtransfoPoly(chipOrder));
}
_mappings[ccdImage->getHashKey()] =
std::make_unique<TwoTransfoMapping>(_chipMap[chip], _visitMap[visit]);
Expand Down

0 comments on commit 430e93e

Please sign in to comment.