30 #include "Eigen/Sparse"
55 double wxx = refStar.
vy /
det;
56 double wyy = refStar.
vx /
det;
57 double wxy = -refStar.
vxy /
det;
58 return wxx *
std::pow(refStar.
x, 2) + 2 * wxy * refStar.
x * refStar.
y + wyy *
std::pow(refStar.
y, 2);
68 _astrometryModel(astrometryModel),
69 _epoch(_associations->getEpoch()),
78 double deltaYears)
const {
83 fittedStarInTP = sky2TP.
apply(temp);
85 fittedStarInTP = sky2TP.
apply(fittedStar);
87 return fittedStarInTP;
93 static void tweakAstromMeasurementErrors(FatPoint &P, MeasuredStar
const &Ms,
double error) {
94 static bool called =
false;
95 static double increment = 0;
106 void AstrometryFit::leastSquareDerivativesMeasurement(CcdImage
const &ccdImage, TripletList &tripletList,
107 Eigen::VectorXd &fullGrad,
108 MeasuredStarList
const *msList)
const {
118 if (msList) assert(&(msList->front()->getCcdImage()) == &ccdImage);
121 const AstrometryMapping *mapping = _astrometryModel->getMapping(ccdImage);
123 std::size_t npar_mapping = (_fittingDistortions) ? mapping->getNpar() : 0;
126 std::size_t npar_tot = npar_mapping + npar_pos + npar_pm;
129 if (npar_tot == 0)
return;
131 if (_fittingDistortions) mapping->getMappingIndices(indices);
134 double deltaYears = _epoch - ccdImage.getEpoch();
136 auto sky2TP = _astrometryModel->getSkyToTangentPlane(ccdImage);
138 AstrometryTransformLinear dypdy;
142 Eigen::Matrix2d transW(2, 2);
143 Eigen::Matrix2d alpha(2, 2);
144 Eigen::VectorXd grad(npar_tot);
146 Eigen::Index kTriplets = tripletList.getNextFreeIndex();
147 const MeasuredStarList &catalog = (msList) ? *msList : ccdImage.getCatalogForFit();
149 for (
auto &i : catalog) {
150 const MeasuredStar &ms = *i;
151 if (!ms.isValid())
continue;
154 tweakAstromMeasurementErrors(inPos, ms, _posError);
158 if (_fittingDistortions)
159 mapping->computeTransformAndDerivatives(inPos, outPos, H);
161 mapping->transformPosAndErrors(inPos, outPos);
164 double det = outPos.vx * outPos.vy -
std::pow(outPos.vxy, 2);
165 if (
det <= 0 || outPos.vx <= 0 || outPos.vy <= 0) {
166 LOGLS_WARN(
_log,
"Inconsistent measurement errors: drop measurement at "
167 << Point(ms) <<
" in image " << ccdImage.getName());
170 transW(0, 0) = outPos.vy /
det;
171 transW(1, 1) = outPos.vx /
det;
172 transW(0, 1) = transW(1, 0) = -outPos.vxy /
det;
175 alpha(0, 0) =
sqrt(transW(0, 0));
177 alpha(1, 0) = transW(0, 1) / alpha(0, 0);
180 alpha(1, 1) = 1. /
sqrt(
det * transW(0, 0));
185 Point fittedStarInTP = transformFittedStar(*fs, *sky2TP, deltaYears);
193 H(npar_mapping, 0) = -dypdy.A11();
194 H(npar_mapping + 1, 0) = -dypdy.A12();
195 H(npar_mapping, 1) = -dypdy.A21();
196 H(npar_mapping + 1, 1) = -dypdy.A22();
197 indices[npar_mapping] = fs->getIndexInMatrix();
198 indices.at(npar_mapping + 1) = fs->getIndexInMatrix() + 1;
214 Eigen::Vector2d res(fittedStarInTP.x - outPos.x, fittedStarInTP.y - outPos.y);
222 for (
std::size_t ipar = 0; ipar < npar_tot; ++ipar) {
224 double val = halpha(ipar, ic);
225 if (
val == 0)
continue;
226 tripletList.addTriplet(indices[ipar], kTriplets + ic,
val);
228 fullGrad(indices[ipar]) += grad(ipar);
232 tripletList.setNextFreeIndex(kTriplets);
235 void AstrometryFit::leastSquareDerivativesReference(FittedStarList
const &fittedStarList,
236 TripletList &tripletList,
237 Eigen::VectorXd &fullGrad)
const {
246 if (!_fittingPos)
return;
250 Eigen::Matrix2d W(2, 2);
251 Eigen::Matrix2d alpha(2, 2);
252 Eigen::Matrix2d H(2, 2), halpha(2, 2), HW(2, 2);
253 AstrometryTransformLinear der;
254 Eigen::Vector2d res, grad;
255 Eigen::Index indices[2];
256 Eigen::Index kTriplets = tripletList.getNextFreeIndex();
263 TanRaDecToPixel proj(AstrometryTransformLinear(), Point(0., 0.));
264 for (
auto const &i : fittedStarList) {
265 const FittedStar &fs = *i;
266 auto rs = fs.getRefStar();
267 if (rs ==
nullptr)
continue;
268 proj.setTangentPoint(fs);
271 proj.transformPosAndErrors(*rs, rsProj);
273 proj.computeDerivative(fs, der, 1e-4);
275 H(0, 0) = -der.A11();
276 H(1, 0) = -der.A12();
277 H(0, 1) = -der.A21();
278 H(1, 1) = -der.A22();
280 double det = rsProj.vx * rsProj.vy -
std::pow(rsProj.vxy, 2);
281 if (rsProj.vx <= 0 || rsProj.vy <= 0 ||
det <= 0) {
282 LOGLS_WARN(
_log,
"RefStar error matrix not positive definite for: " << *rs);
285 W(0, 0) = rsProj.vy /
det;
286 W(0, 1) = W(1, 0) = -rsProj.vxy /
det;
287 W(1, 1) = rsProj.vx /
det;
290 alpha(0, 0) =
sqrt(W(0, 0));
292 alpha(1, 0) = W(0, 1) / alpha(0, 0);
293 alpha(1, 1) = 1. /
sqrt(
det * W(0, 0));
295 indices[0] = fs.getIndexInMatrix();
296 indices[1] = fs.getIndexInMatrix() + 1;
297 unsigned npar_tot = 2;
313 for (
std::size_t ipar = 0; ipar < npar_tot; ++ipar) {
314 for (
unsigned ic = 0; ic < 2; ++ic) {
315 double val = halpha(ipar, ic);
316 if (
val == 0)
continue;
317 tripletList.addTriplet(indices[ipar], kTriplets + ic,
val);
319 fullGrad(indices[ipar]) += grad(ipar);
323 tripletList.setNextFreeIndex(kTriplets);
326 void AstrometryFit::accumulateStatImage(CcdImage
const &ccdImage, Chi2Accumulator &accum)
const {
333 const AstrometryMapping *mapping = _astrometryModel->getMapping(ccdImage);
335 double deltaYears = _epoch - ccdImage.getEpoch();
337 auto sky2TP = _astrometryModel->getSkyToTangentPlane(ccdImage);
339 Eigen::Matrix2Xd transW(2, 2);
341 auto &catalog = ccdImage.getCatalogForFit();
342 for (
auto const &ms : catalog) {
343 if (!ms->isValid())
continue;
345 FatPoint inPos = *ms;
346 tweakAstromMeasurementErrors(inPos, *ms, _posError);
350 mapping->transformPosAndErrors(inPos, outPos);
351 double det = outPos.vx * outPos.vy -
std::pow(outPos.vxy, 2);
352 if (
det <= 0 || outPos.vx <= 0 || outPos.vy <= 0) {
353 LOGLS_WARN(
_log,
" Inconsistent measurement errors :drop measurement at "
354 << Point(*ms) <<
" in image " << ccdImage.getName());
357 transW(0, 0) = outPos.vy /
det;
358 transW(1, 1) = outPos.vx /
det;
359 transW(0, 1) = transW(1, 0) = -outPos.vxy /
det;
362 Point fittedStarInTP = transformFittedStar(*fs, *sky2TP, deltaYears);
364 Eigen::Vector2d res(fittedStarInTP.x - outPos.x, fittedStarInTP.y - outPos.y);
365 double chi2Val = res.transpose() * transW * res;
367 accum.addEntry(chi2Val, 2, ms);
371 void AstrometryFit::accumulateStatImageList(
CcdImageList const &ccdImageList, Chi2Accumulator &accum)
const {
372 for (
auto const &ccdImage : ccdImageList) {
373 accumulateStatImage(*ccdImage, accum);
377 void AstrometryFit::accumulateStatRefStars(Chi2Accumulator &accum)
const {
385 FittedStarList &fittedStarList =
_associations->fittedStarList;
386 TanRaDecToPixel proj(AstrometryTransformLinear(), Point(0., 0.));
387 for (
auto const &fs : fittedStarList) {
388 auto rs = fs->getRefStar();
389 if (rs ==
nullptr)
continue;
390 proj.setTangentPoint(*fs);
393 proj.transformPosAndErrors(*rs, rsProj);
395 double chi2 = computeProjectedRefStarChi2(rsProj);
396 accum.addEntry(chi2, 2, fs);
403 void AstrometryFit::getIndicesOfMeasuredStar(MeasuredStar
const &measuredStar,
IndexVector &indices)
const {
404 if (_fittingDistortions) {
405 const AstrometryMapping *mapping = _astrometryModel->getMapping(measuredStar.getCcdImage());
406 mapping->getMappingIndices(indices);
409 Eigen::Index fsIndex = fs->getIndexInMatrix();
425 _fittingDistortions = (
_whatToFit.
find(
"Distortions") != std::string::npos);
426 _fittingPos = (
_whatToFit.
find(
"Positions") != std::string::npos);
436 for (
auto &fittedStar : fittedStarList) {
456 "AstrometryFit::offsetParams : the provided vector length is not compatible with "
457 "the current whatToFit setting");
458 if (_fittingDistortions) _astrometryModel->offsetParams(delta);
462 for (
auto const &i : fittedStarList) {
468 fs.
x += delta(index);
469 fs.
y += delta(index + 1);
480 const char *what2fit[] = {
"Positions",
"Distortions",
"Positions Distortions"};
482 for (
unsigned k = 0; k <
sizeof(what2fit) /
sizeof(what2fit[0]); ++k) {
489 jacobian.setFromTriplets(tripletList.
begin(), tripletList.
end());
499 ofile <<
"id" << separator <<
"xccd" << separator <<
"yccd " << separator;
500 ofile <<
"rx" << separator <<
"ry" << separator;
501 ofile <<
"xtp" << separator <<
"ytp" << separator;
502 ofile <<
"mag" << separator <<
"deltaYears" << separator;
503 ofile <<
"xErr" << separator <<
"yErr" << separator <<
"xyCov" << separator;
504 ofile <<
"xtpi" << separator <<
"ytpi" << separator;
505 ofile <<
"rxi" << separator <<
"ryi" << separator;
506 ofile <<
"fsindex" << separator;
507 ofile <<
"ra" << separator <<
"dec" << separator;
508 ofile <<
"chi2" << separator <<
"nm" << separator;
509 ofile <<
"detector" << separator <<
"visit" <<
std::endl;
511 ofile <<
"id in source catalog" << separator <<
"coordinates in CCD (pixels)" << separator << separator;
512 ofile <<
"residual on TP (degrees)" << separator << separator;
513 ofile <<
"transformed coordinate in TP (degrees)" << separator << separator;
514 ofile <<
"rough magnitude" << separator <<
"Julian epoch year delta from fit epoch" << separator;
515 ofile <<
"transformed measurement uncertainty (degrees)" << separator << separator << separator;
516 ofile <<
"as-read position on TP (degrees)" << separator << separator;
517 ofile <<
"as-read residual on TP (degrees)" << separator << separator;
518 ofile <<
"unique index of the fittedStar" << separator;
519 ofile <<
"on sky position of fittedStar" << separator << separator;
520 ofile <<
"contribution to chi2 (2D dofs)" << separator <<
"number of measurements of this fittedStar"
522 ofile <<
"detector id" << separator <<
"visit id" <<
std::endl;
525 for (
auto const &ccdImage : ccdImageList) {
528 const auto readTransform = ccdImage->getReadWcs();
530 double deltaYears = _epoch - ccdImage->getEpoch();
531 for (
auto const &ms : cat) {
532 if (!ms->isValid())
continue;
535 tweakAstromMeasurementErrors(inPos, *ms, _posError);
537 auto sky2TP = _astrometryModel->getSkyToTangentPlane(*ccdImage);
539 compose(*sky2TP, *readTransform);
540 FatPoint inputTpPos = readPixToTangentPlane->apply(inPos);
543 Point fittedStarInTP = transformFittedStar(*fs, *sky2TP, deltaYears);
544 Point res = tpPos - fittedStarInTP;
545 Point inputRes = inputTpPos - fittedStarInTP;
547 double wxx = tpPos.
vy /
det;
548 double wyy = tpPos.
vx /
det;
549 double wxy = -tpPos.
vxy /
det;
550 double chi2 = wxx * res.
x * res.
x + wyy * res.
y * res.
y + 2 * wxy * res.
x * res.
y;
553 ofile << ms->getId() << separator << ms->x << separator << ms->y << separator;
554 ofile << res.
x << separator << res.
y << separator;
555 ofile << tpPos.
x << separator << tpPos.
y << separator;
556 ofile << fs->getMag() << separator << deltaYears << separator;
557 ofile << tpPos.
vx << separator << tpPos.
vy << separator << tpPos.
vxy << separator;
558 ofile << inputTpPos.
x << separator << inputTpPos.
y << separator;
559 ofile << inputRes.
x << separator << inputRes.
y << separator;
560 ofile << fs->getIndexInMatrix() << separator;
561 ofile << fs->x << separator << fs->y << separator;
562 ofile << chi2 << separator << fs->getMeasurementCount() << separator;
563 ofile << ccdImage->getCcdId() << separator << ccdImage->getVisit() <<
std::endl;
572 ofile <<
"ra" << separator <<
"dec " << separator;
573 ofile <<
"rx" << separator <<
"ry" << separator;
574 ofile <<
"mag" << separator;
575 ofile <<
"xErr" << separator <<
"yErr" << separator <<
"xyCov" << separator;
576 ofile <<
"fsindex" << separator;
577 ofile <<
"chi2" << separator <<
"nm" <<
std::endl;
579 ofile <<
"coordinates of fittedStar (degrees)" << separator << separator;
580 ofile <<
"residual on TP (degrees)" << separator << separator;
581 ofile <<
"magnitude" << separator;
582 ofile <<
"refStar transformed measurement uncertainty (degrees)" << separator << separator << separator;
583 ofile <<
"unique index of the fittedStar" << separator;
584 ofile <<
"refStar contribution to chi2 (2D dofs)" << separator
585 <<
"number of measurements of this FittedStar" <<
std::endl;
590 for (
auto const &i : fittedStarList) {
593 if (rs ==
nullptr)
continue;
598 double chi2 = computeProjectedRefStarChi2(rsProj);
601 ofile << fs.
x << separator << fs.
y << separator;
602 ofile << rsProj.
x << separator << rsProj.
y << separator;
603 ofile << fs.
getMag() << separator;
604 ofile << rsProj.
vx << separator << rsProj.
vy << separator << rsProj.
vxy << separator;
Eigen::Matrix< double, Eigen::Dynamic, 2 > MatrixX2d
Eigen::SparseMatrix< double, 0, Eigen::Index > SparseMatrixD
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
LSST DM logging module built on log4cxx.
#define LOGLS_WARN(logger, message)
Log a warn-level message using an iostream-based interface.
#define LOG_GET(logger)
Returns a Log object associated with logger.
#define LOGLS_INFO(logger, message)
Log a info-level message using an iostream-based interface.
#define LOGLS_DEBUG(logger, message)
Log a debug-level message using an iostream-based interface.
void checkStuff()
DEBUGGING routine.
void saveChi2RefContributions(std::string const &filename) const override
Save a CSV file containing residuals of reference terms.
void offsetParams(Eigen::VectorXd const &delta) override
Offset the parameters by the requested quantities.
AstrometryFit(std::shared_ptr< Associations > associations, std::shared_ptr< AstrometryModel > astrometryModel, double posError)
this is the only constructor
void assignIndices(std::string const &whatToFit) override
Set parameters to fit and assign indices in the big matrix.
void saveChi2MeasContributions(std::string const &filename) const override
Save a CSV file containing residuals of measurement terms.
virtual class needed in the abstraction of the distortion model
virtual void transformPosAndErrors(FatPoint const &where, FatPoint &outPoint) const =0
The same as above but without the parameter derivatives (used to evaluate chi^2)
A Point with uncertainties.
FittedStars are objects whose position or flux is going to be fitted, and which come from the associa...
const RefStar * getRefStar() const
Get the astrometric reference star associated with this star.
void setIndexInMatrix(Eigen::Index const index)
index is a value that a fit can set and reread....
int getMeasurementCount() const
The number of MeasuredStars currently associated with this FittedStar.
Eigen::Index getIndexInMatrix() const
A list of FittedStar s. Such a list is typically constructed by Associations.
void leastSquareDerivatives(TripletList &tripletList, Eigen::VectorXd &grad) const
Evaluates the chI^2 derivatives (Jacobian and gradient) for the current whatToFit setting.
Eigen::Index _nModelParams
std::shared_ptr< Associations > _associations
Eigen::Index _nStarParams
A list of MeasuredStar. They are usually filled in Associations::createCcdImage.
Point applyProperMotion(Point star, double timeDeltaYears) const
Apply proper motion correction to the input star, returning a star with PM-corrected coordinates and ...
This one is the Tangent Plane (called gnomonic) projection (from celestial sphere to tangent plane)
void transformPosAndErrors(const FatPoint &in, FatPoint &out) const
transform with analytical derivatives
void setTangentPoint(Point const &tangentPoint)
Resets the projection (or tangent) point.
Eigen::Index getNextFreeIndex() const
Reports invalid arguments.
std::list< std::shared_ptr< CcdImage > > CcdImageList
std::unique_ptr< AstrometryTransform > compose(AstrometryTransform const &left, AstrometryTransform const &right)
Returns a pointer to a composition of transforms, representing left(right()).
A base class for image defects.
T setprecision(T... args)