78 centers.
push_back(ccdImage->getBoresightRaDec());
86 _commonTangentPoint =
Point(commonTangentPoint.getX(), commonTangentPoint.getY());
87 for (
auto &ccdImage :
ccdImageList) ccdImage->setCommonTangentPoint(_commonTangentPoint);
92 Frame tangentPlaneFrame;
95 Frame CTPFrame = ccdImage->getPixelToCommonTangentPlane()->apply(ccdImage->getImageFrame(),
false);
96 if (tangentPlaneFrame.
getArea() == 0)
97 tangentPlaneFrame = CTPFrame;
99 tangentPlaneFrame += CTPFrame;
104 TanPixelToRaDec commonTangentPlaneToRaDec(identity, _commonTangentPoint);
105 Frame raDecFrame = commonTangentPlaneToRaDec.
apply(tangentPlaneFrame,
false);
120 const bool enlargeFittedList) {
126 item->clearBeforeAssoc();
136 ccdImage->resetCatalogForFit();
142 Frame ccdImageFrameCPT = toCommonTangentPlane->apply(ccdImage->getImageFrame(),
false);
143 ccdImageFrameCPT = ccdImageFrameCPT.
rescale(1.10);
149 if (ccdImageFrameCPT.
inFrame(*fittedStar)) {
156 toCommonTangentPlane.
get(), matchCutInArcSec / 3600.);
159 LOGLS_DEBUG(_log,
"Measured-to-Fitted matches before removing ambiguities " << starMatchList->size());
160 starMatchList->removeAmbiguities(*toCommonTangentPlane);
161 LOGLS_DEBUG(_log,
"Measured-to-Fitted matches after removing ambiguities " << starMatchList->size());
165 int matchedCount = 0;
166 for (
auto const &starMatch : *starMatchList) {
167 auto bs = starMatch.s1;
168 auto ms_const = std::dynamic_pointer_cast<const MeasuredStar>(bs);
169 auto ms = std::const_pointer_cast<MeasuredStar>(ms_const);
170 auto bs2 = starMatch.s2;
171 auto fs_const = std::dynamic_pointer_cast<const FittedStar>(bs2);
172 auto fs = std::const_pointer_cast<FittedStar>(fs_const);
173 ms->setFittedStar(fs);
176 LOGLS_DEBUG(_log,
"Matched " << matchedCount <<
" objects in " << ccdImage->getName());
179 int unMatchedCount = 0;
180 for (
auto const &mstar : catalog) {
182 if (mstar->getFittedStar())
continue;
183 if (enlargeFittedList) {
184 auto fs = std::make_shared<FittedStar>(*mstar);
186 toCommonTangentPlane->transformPosAndErrors(*fs, *fs);
188 mstar->setFittedStar(fs);
192 LOGLS_DEBUG(_log,
"Unmatched objects: " << unMatchedCount);
204 std::string const &fluxField,
float refCoordinateErr,
205 bool rejectBadFluxes) {
206 if (refCat.
size() == 0) {
208 " reference catalog is empty : stop here "));
216 raErrKey = refCat.
getSchema()[
"coord_raErr"];
217 decErrKey = refCat.
getSchema()[
"coord_decErr"];
220 auto fluxKey = refCat.
getSchema().
find<
double>(fluxField).key;
224 fluxErrKey = refCat.
getSchema().
find<
double>(fluxField +
"Err").key;
227 << fluxField <<
"Err"
228 <<
") not found in reference catalog. Not using ref flux errors.");
238 LOGLS_WARN(_log,
"Not loading proper motions: (pm_ra,pm_dec) fields not found in reference catalog.");
240 LOGLS_WARN(_log,
"Not loading proper motions: RA/Dec proper motion values must be `geom:Angle`: "
245 pmDecErrKey = refCat.
getSchema().
find<
float>(
"pm_decErr").key;
247 LOGLS_WARN(_log,
"Not loading proper motions: error fields not available: " <<
ex.what());
252 pmRaDecCovKey = refCat.
getSchema().
find<
float>(
"pm_ra_Dec_Cov").key;
254 LOGLS_WARN(_log,
"No ra/dec proper motion covariances in refcat: " <<
ex.what());
258 for (
size_t i = 0; i < refCat.
size(); i++) {
259 auto const &record = refCat.
get(i);
261 auto coord = record->
get(coordKey);
262 double flux = record->get(fluxKey);
264 if (fluxErrKey.isValid()) {
265 fluxErr = record->get(fluxErrKey);
271 auto star = std::make_shared<RefStar>(ra,
dec, flux, fluxErr);
279 star->vx =
std::pow(refCoordinateErr / 1000. / 3600. /
std::cos(coord.getLatitude()), 2);
280 star->vy =
std::pow(refCoordinateErr / 1000. / 3600., 2);
285 star->setProperMotion(std::make_unique<ProperMotion const>(
286 record->get(pmRaKey).asRadians(), record->get(pmDecKey).asRadians(),
287 record->get(pmRaErrKey), record->get(pmDecErrKey), record->get(pmRaDecCovKey)));
289 star->setProperMotion(std::make_unique<ProperMotion const>(
290 record->get(pmRaKey).asRadians(), record->get(pmDecKey).asRadians(),
291 record->get(pmRaErrKey), record->get(pmDecErrKey)));
305 TanRaDecToPixel raDecToCommonTangentPlane(identity, _commonTangentPoint);
307 associateRefStars(matchCut.
asArcseconds(), &raDecToCommonTangentPlane);
314 matchCutInArcSec / 3600.);
316 LOGLS_DEBUG(_log,
"Refcat matches before removing ambiguities " << starMatchList->size());
317 starMatchList->removeAmbiguities(*
transform);
318 LOGLS_DEBUG(_log,
"Refcat matches after removing ambiguities " << starMatchList->size());
321 for (
auto const &starMatch : *starMatchList) {
325 const BaseStar &bs2 = *starMatch.s2;
333 "Associated " << starMatchList->size() <<
" reference stars among " <<
refStarList.
size());
337 selectFittedStars(minMeasurements);
338 normalizeFittedStars();
341 void Associations::selectFittedStars(
int minMeasurements) {
344 int totalMeasured = 0, validMeasured = 0;
356 if (fittedStar ==
nullptr) {
363 if (!fittedStar->getRefStar() && fittedStar->getMeasurementCount() < minMeasurements) {
364 fittedStar->getMeasurementCount()--;
365 mi = catalog.
erase(mi);
375 if ((*fi)->getMeasurementCount() == 0) {
383 LOGLS_INFO(_log,
"Total, valid number of Measured stars: " << totalMeasured <<
", " << validMeasured);
386 void Associations::normalizeFittedStars() {
391 fittedStar->setFlux(0.0);
392 fittedStar->getMag() = 0.0;
398 MeasuredStarList &catalog = ccdImage->getCatalogForFit();
399 for (
auto &mi : catalog) {
400 auto fittedStar = mi->getFittedStar();
401 if (fittedStar ==
nullptr)
403 pex::exceptions::RuntimeError,
404 "All measuredStars must have a fittedStar: did you call selectFittedStars()?"));
405 auto point = toCommonTangentPlane->apply(*mi);
406 fittedStar->x += point.x;
407 fittedStar->y += point.y;
408 fittedStar->getFlux() += mi->getFlux();
412 _maxMeasuredStars = 0;
414 auto measurementCount = fi->getMeasurementCount();
415 _maxMeasuredStars += measurementCount;
416 fi->x /= measurementCount;
417 fi->y /= measurementCount;
418 fi->getFlux() /= measurementCount;
428 auto fittedStar = *
iter;
429 if (fittedStar->getMeasurementCount() == 0) {
430 LOGLS_TRACE(_log,
"Deleting FittedStar (has no measuredStars): " << *fittedStar);
438 LOGLS_INFO(_log,
"Removed " <<
count <<
" fittedStars that had no associated measuredStar.");
442 void Associations::assignMags() {
445 for (
auto const &mstar : catalog) {
447 if (!fstar)
continue;
458 "DeprojectFittedStars: Fitted stars are already in sidereal coordinates, nothing done ");
469 return item->getCatalogForFit().size() > 0;
476 if ((fittedStar !=
nullptr) & (fittedStar->getRefStar() !=
nullptr))
count++;
482 void Associations::collectMCStars(
int realization) {
488 string dbimdir = ccdImage.Dir();
489 string mctruth = dbimdir +
"/mc/mctruth.list";
491 if (realization >= 0) {
493 sstrm << dbimdir <<
"/mc/mctruth_" << realization <<
".list";
494 mctruth = sstrm.
str();
497 AstrometryTransformIdentity gti;
501 DicStarList mctruthlist(mctruth);
505 for (smI = starMatchList->begin(); smI != starMatchList->end(); smI++) {
506 StarMatch &sm = *smI;
507 BaseStar *bs = sm.s1;
508 MeasuredStar *mstar =
dynamic_cast<MeasuredStar *
>(bs);
510 DicStar *dstar =
dynamic_cast<DicStar *
>(bs);
512 mcstar->GetMCInfo().iflux = dstar->getval(
"iflux");
513 mcstar->GetMCInfo().tflux = dstar->getval(
"sflux");
520 LOGLS_FATAL(_log,
"CollectMCStars Unable to match MCTruth w/ catalog!");
This file declares a class for representing circular regions on the unit sphere.
This file declares a class for representing convex polygons with great circle edges on the unit spher...
table::Key< int > detector
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Combinatorial searches for linear transformations to go from list1 to list2.
LSST DM logging module built on log4cxx.
#define LOGLS_WARN(logger, message)
Log a warn-level message using an iostream-based interface.
#define LOGLS_FATAL(logger, message)
Log a fatal-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.
#define LOGLS_TRACE(logger, message)
Log a trace-level message using an iostream-based interface.
This file contains a class representing spherical coordinates.
table::Key< table::Array< std::uint8_t > > wcs
size_type size() const
Return the number of elements in the catalog.
Schema getSchema() const
Return the schema associated with the catalog's table.
std::shared_ptr< RecordT > const get(size_type i) const
Return a pointer to the record at index i.
A FunctorKey used to get or set celestial coordinates from a pair of lsst::geom::Angle keys.
bool isValid() const noexcept
Return true if the key was initialized to valid offset.
SchemaItem< T > find(std::string const &name) const
Find a SchemaItem in the Schema by name.
Custom catalog class for record/table subclasses that are guaranteed to have an ID,...
A class representing an angle.
constexpr double asArcseconds() const noexcept
Return an Angle's value in arcseconds.
An integer coordinate rectangle.
size_t nFittedStarsWithAssociatedRefStar() const
Return the number of fittedStars that have an associated refStar.
void cleanFittedStars()
Remove FittedStars that have no measured stars; this can happen after outlier rejection.
void computeCommonTangentPoint()
Sets a shared tangent point for all ccdImages, using the mean of the centers of all ccdImages.
lsst::sphgeom::Circle computeBoundingCircle() const
Return the bounding circle in on-sky (RA, Dec) coordinates containing all CcdImages.
int nCcdImagesValidForFit() const
return the number of CcdImages with non-empty catalogs to-be-fit.
void createCcdImage(afw::table::SourceCatalog &catalog, std::shared_ptr< lsst::afw::geom::SkyWcs > wcs, std::shared_ptr< lsst::afw::image::VisitInfo > visitInfo, lsst::geom::Box2I const &bbox, std::string const &filter, std::shared_ptr< afw::image::PhotoCalib > photoCalib, std::shared_ptr< afw::cameraGeom::Detector > detector, int visit, int ccd, lsst::jointcal::JointcalControl const &control)
Create a ccdImage from an exposure catalog and metadata, and add it to the list.
void setCommonTangentPoint(lsst::geom::Point2D const &commonTangentPoint)
Shared tangent point for all ccdImages (decimal degrees).
void associateCatalogs(const double matchCutInArcsec=0, const bool useFittedList=false, const bool enlargeFittedList=true)
incrementaly builds a merged catalog of all image catalogs
Point getCommonTangentPoint() const
Shared tangent point for all ccdImages (decimal degrees).
void collectRefStars(afw::table::SimpleCatalog &refCat, geom::Angle matchCut, std::string const &fluxField, float refCoordinateErr, bool rejectBadFluxes=false)
Collect stars from an external reference catalog and associate them with fittedStars.
CcdImageList ccdImageList
FittedStarList fittedStarList
void prepareFittedStars(int minMeasurements)
Prepare the fittedStar list by making quality cuts and normalizing measurements.
void deprojectFittedStars()
Sends back the fitted stars coordinates on the sky FittedStarsList::inTangentPlaneCoordinates keeps t...
The base class for handling stars. Used by all matching routines.
Handler of an actual image from a single CCD.
MeasuredStarList const & getCatalogForFit() const
Gets the catalog to be used for fitting, which may have been cleaned-up.
FittedStars are objects whose position or flux is going to be fitted, and which come from the associa...
void setRefStar(const RefStar *_refStar)
Set the astrometric reference star associated with this star.
A list of FittedStar s. Such a list is typically constructed by Associations.
bool inTangentPlaneCoordinates
rectangle with sides parallel to axes.
Frame rescale(const double factor) const
rescale it. The center does not move.
bool inFrame(double x, double y) const
inside?
double xMin
coordinate of boundary.
Sources measured on images.
std::shared_ptr< FittedStar > getFittedStar() const
double getMagWeight() const
the inverse of the mag variance
A list of MeasuredStar. They are usually filled in Associations::createCcdImage.
Objects used as position/flux anchors (e.g.
void applyTransform(const Operator &op)
enables to apply a geometrical transform if Star is Basestar or derives from it.
The transformation that handles pixels to sideral transformations (Gnomonic, possibly with polynomial...
virtual void apply(const double xIn, const double yIn, double &xOut, double &yOut) const=0
This one is the Tangent Plane (called gnomonic) projection (from celestial sphere to tangent plane)
Reports invalid arguments.
Reports attempts to access elements using an invalid key.
Reports errors from accepting an object of an unexpected or inappropriate type.
Circle is a circular region on the unit sphere that contains its boundary.
static ConvexPolygon convexHull(std::vector< UnitVector3d > const &points)
convexHull returns the convex hull of the given set of points if it exists and throws an exception ot...
Circle getBoundingCircle() const override
getBoundingCircle returns a bounding-circle for this region.
static LonLat fromDegrees(double lon, double lat)
T emplace_back(T... args)
double nanojanskyToABMagnitude(double flux)
Convert a flux in nanojansky to AB magnitude.
SpherePoint averageSpherePoint(std::vector< SpherePoint > const &coords)
Return the average of a list of coordinates.
constexpr AngleUnit degrees
constant with units of degrees
constexpr double radToDeg(double x) noexcept
BaseStarList & Measured2Base(MeasuredStarList &This)
std::unique_ptr< StarMatchList > listMatchCollect(const BaseStarList &list1, const BaseStarList &list2, const AstrometryTransform *guess, const double maxDist)
assembles star matches.
BaseStarList & Fitted2Base(FittedStarList &This)
MeasuredStarList::iterator MeasuredStarIterator
::std::list< StarMatch >::iterator StarMatchIterator
BaseStarList & Ref2Base(RefStarList &This)
FittedStarList::iterator FittedStarIterator
A base class for image defects.
std::string sourceFluxField
"name of flux field in source catalog" ;