Skip to content

Commit

Permalink
Compute common tangent plane after all ccdImages are loaded
Browse files Browse the repository at this point in the history
  • Loading branch information
@parejkoj
parejkoj committed Feb 10, 2017
1 parent 407a09a commit b0ecb53
Show file tree
Hide file tree
Showing 5 changed files with 91 additions and 70 deletions.
34 changes: 23 additions & 11 deletions include/lsst/jointcal/Associations.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,8 @@ class Associations {
RefStarList refStarList; // e.g. GAIA or SDSS reference stars
FittedStarList fittedStarList; // stars that are going to be fitted

Point _commonTangentPoint;
public:
// TODO: NOTE: these only exist because those lists can't be swig'd because
// swig doesn't like boost::intrusive_ptr (which the lists contain).
// Once DM-4043 is solved, delete these.
Expand All @@ -42,28 +44,37 @@ class Associations {
size_t fittedStarListSize()
{ return fittedStarList.size(); }

// fit cuts and stuff:
Point _commonTangentPoint;

void assignMags();

public:

/**
* Source selection is performed in python, so Associations' constructor
* only initializes a couple of variables.
*/
Associations();

//! Sets a tangent point (reasonably centered for the input image set).
void setCommonTangentPoint(Point const &commonTangentPoint)
{ _commonTangentPoint = commonTangentPoint;};
/**
* @brief Sets a shared tangent point for all ccdImages.
*
* @param commonTangentPoint The common tangent point of all input images (decimal degrees).
*/
void setCommonTangentPoint(lsst::afw::geom::Point2D const &commonTangentPoint);

//! can be used to project sidereal coordinates related to the image set on a plane.
Point getCommonTangentPoint() const { return _commonTangentPoint;}

//! same as above
bool addImage(lsst::afw::table::SortedCatalogT<lsst::afw::table::SourceRecord> &catalog,
/**
* @brief Create a ccdImage from an exposure catalog and metadata, and add it to the list.
*
* @param catalog The extracted source catalog, selected for good astrometric sources.
* @param[in] wcs The exposure's original wcs
* @param[in] visitInfo The exposure's visitInfo object
* @param bbox The bounding box of the exposure
* @param filter The exposure's filter
* @param[in] calib The exposure's photometric calibration
* @param[in] visit The visit identifier
* @param[in] ccd The ccd identifier
* @param[in] control The JointcalControl object
*/
void addImage(lsst::afw::table::SortedCatalogT<lsst::afw::table::SourceRecord> &catalog,
std::shared_ptr<lsst::afw::image::TanWcs> wcs,
std::shared_ptr<lsst::afw::image::VisitInfo> visitInfo,
lsst::afw::geom::Box2I const &bbox,
Expand Down Expand Up @@ -109,6 +120,7 @@ class Associations {
private:
void associateRefStars(double matchCutInArcsec, const Gtransfo* gtransfo);

void assignMags();
};

#endif /* ASSOCIATIONS__H */
Expand Down
29 changes: 22 additions & 7 deletions include/lsst/jointcal/CcdImage.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,13 +9,12 @@
#include "lsst/afw/image/TanWcs.h"
#include "lsst/afw/image/Calib.h"
#include "lsst/afw/image/VisitInfo.h"
#include "lsst/afw/coord/Coord.h"
#include "lsst/daf/base/PropertySet.h"
#include "lsst/afw/geom/Box.h"
#include "lsst/jointcal/MeasuredStar.h"
#include "lsst/jointcal/Gtransfo.h"
#include "lsst/jointcal/Frame.h"
//#include "lsst/jointcal/Jointcal.h"


namespace lsst {
namespace jointcal {
Expand Down Expand Up @@ -51,6 +50,7 @@ class CcdImage : public RefCount
CcdIdType _ccdId;
VisitIdType _visit;

lsst::afw::coord::IcrsCoord boresightRaDec;
double airMass; // airmass value.
double fluxCoeff; // coefficient to convert ADUs to ADUs/sec at airmass 1
double mjd; // modified julian date
Expand All @@ -61,14 +61,13 @@ class CcdImage : public RefCount

std::string _filter;

Point commonTangentPoint;
Point _commonTangentPoint;

void LoadCatalog(const lsst::afw::table::SortedCatalogT<lsst::afw::table::SourceRecord> &Cat, const std::string &fluxField);

public:

CcdImage(lsst::afw::table::SortedCatalogT<lsst::afw::table::SourceRecord> &Ri,
const Point &CommonTangentPoint,
const PTR(lsst::afw::image::TanWcs) wcs,
const PTR(lsst::afw::image::VisitInfo) visitInfo,
const lsst::afw::geom::Box2I &bbox,
Expand Down Expand Up @@ -98,6 +97,20 @@ class CcdImage : public RefCount
MeasuredStarList & getCatalogForFit() { return _catalogForFit;}
//@}

/**
* @brief Sets the common tangent point and computes necessary transforms.
*
* @param[in] commonTangentPoint The common tangent point of all ccdImages (decimal degrees).
*/
void setCommonTangentPoint(const Point &commonTangentPoint);

/**
* @brief Gets the common tangent point, shared between all ccdImages.
*
* @return The common tangent point of all ccdImages (decimal degrees).
*/
Point const& getCommonTangentPoint() const { return _commonTangentPoint; }

//!
const Gtransfo* Pix2CommonTangentPlane() const
{ return pix2CommonTangentPlane.get();}
Expand Down Expand Up @@ -136,6 +149,11 @@ class CcdImage : public RefCount
//!Return the exposure's photometric calibration
PTR(lsst::afw::image::Calib) getCalib() { return _calib; }

/**
* @brief Gets the boresight RA/Dec.
*/
lsst::afw::coord::IcrsCoord getBoresightRaDec() { return boresightRaDec; }

//!
double HourAngle() const { return hourAngle; }

Expand Down Expand Up @@ -166,9 +184,6 @@ class CcdImage : public RefCount
//! Frame in pixels
const Frame& ImageFrame() const { return imageFrame;}

//! Common Tangent Point
Point const& CommonTangentPoint() const { return commonTangentPoint; }

private:
CcdImage(const CcdImage &); // forbid copies
};
Expand Down
5 changes: 5 additions & 0 deletions python/lsst/jointcal/jointcal.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -259,6 +259,11 @@ def run(self, dataRefs, profile_jointcal=False):
filters.append(result.filter)
filters = collections.Counter(filters)

centers = [ccdImage.getBoresightRaDec() for ccdImage in associations.getCcdImageList()]
commonTangentPoint = lsst.afw.coord.averageCoord(centers)
print("Using common tangent point: ", commonTangentPoint.getPosition())
associations.setCommonTangentPoint(commonTangentPoint.getPosition())

# Use external reference catalogs handled by LSST stack mechanism
# Get the bounding box overlapping all associated images
# ==> This is probably a bad idea to do it this way <== To be improved
Expand Down
27 changes: 9 additions & 18 deletions src/Associations.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ Associations::Associations()
_commonTangentPoint = Point(0, 0);
}

bool Associations::addImage(lsst::afw::table::SortedCatalogT<lsst::afw::table::SourceRecord> &catalog,
void Associations::addImage(lsst::afw::table::SortedCatalogT<lsst::afw::table::SourceRecord> &catalog,
std::shared_ptr<lsst::afw::image::TanWcs> wcs,
std::shared_ptr<lsst::afw::image::VisitInfo> visitInfo,
lsst::afw::geom::Box2I const &bbox,
Expand All @@ -48,26 +48,17 @@ bool Associations::addImage(lsst::afw::table::SortedCatalogT<lsst::afw::table::S
int ccd,
std::shared_ptr<lsst::jointcal::JointcalControl> control)
{
// TODO: I don't like the commonTangentPoint stuff here:
// 1. should create ccdImage first, and extract ra/dec from it
// 2. why take the common point from the first image only? Shouldn't we update it as new images come in?

/* if commonTangentPoint was never initialized
take the one from this image */
if (_commonTangentPoint.x == 0. && _commonTangentPoint.y == 0.)
{
lsst::daf::base::PropertyList::Ptr wcsMeta = wcs->getFitsMetadata();
double crval1 = wcsMeta->get<double>("CRVAL1");
double crval2 = wcsMeta->get<double>("CRVAL2");
_commonTangentPoint = Point(crval1, crval2);
std::cout << "setting common TangentPoint" << _commonTangentPoint << std::endl;
}

std::shared_ptr<CcdImage> ccdImage(new CcdImage(catalog, _commonTangentPoint, wcs, visitInfo, bbox, filter, calib, visit, ccd, control->sourceFluxField));
std::shared_ptr<CcdImage> ccdImage(new CcdImage(catalog, wcs, visitInfo, bbox, filter, calib, visit, ccd, control->sourceFluxField));
ccdImageList.push_back(ccdImage);
std::cout << " we have " << ccdImage->getWholeCatalog().size()
<< " objects in this catalog " << visit << " " << ccd << std::endl;
return true;
}

void Associations::setCommonTangentPoint(lsst::afw::geom::Point2D const &commonTangentPoint)
{
_commonTangentPoint = Point(commonTangentPoint.getX(), commonTangentPoint.getY()); // a jointcal::Point
for (auto &ccdImage: ccdImageList)
ccdImage->setCommonTangentPoint(_commonTangentPoint);
}

void Associations::associateCatalogs(const double matchCutInArcSec,
Expand Down
66 changes: 32 additions & 34 deletions src/CcdImage.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,6 @@ void CcdImage::LoadCatalog(const lsst::afw::table::SortedCatalogT<lsst::afw::tab
}

CcdImage::CcdImage(lsst::afw::table::SortedCatalogT<lsst::afw::table::SourceRecord> &Ri,
const Point &CommonTangentPoint,
const PTR(lsst::afw::image::TanWcs) wcs,
const PTR(lsst::afw::image::VisitInfo) visitInfo,
const lsst::afw::geom::Box2I &bbox,
Expand All @@ -73,8 +72,7 @@ CcdImage::CcdImage(lsst::afw::table::SortedCatalogT<lsst::afw::table::SourceReco
const int &ccdId,
const std::string &fluxField ) :

_filter(filter), _visit(visit), _ccdId(ccdId),
commonTangentPoint(CommonTangentPoint), _calib(calib)
_filter(filter), _visit(visit), _ccdId(ccdId), _calib(calib)

{
LoadCatalog(Ri, fluxField);
Expand All @@ -84,45 +82,18 @@ CcdImage::CcdImage(lsst::afw::table::SortedCatalogT<lsst::afw::table::SourceReco
imageFrame = Frame(lowerLeft, upperRight);

readWcs = new jointcal::TanSipPix2RaDec(jointcal::ConvertTanWcs(wcs));

// use some other variable in case we later have to actually convert the
// read wcs:
const BaseTanWcs* tanWcs = readWcs.get();

inverseReadWcs = readWcs->InverseTransfo(0.01, imageFrame);

std::stringstream out;
out << visit << "_" << ccdId;
name = out.str();

/* we don't assume here that we know the internals of TanPix2RaDec:
to construct pix->TP, we do pix->sky->TP, although pix->sky
actually goes through TP */

GtransfoLin identity;
TanRaDec2Pix raDec2TP(identity, tanWcs->TangentPoint());
pix2TP = GtransfoCompose(&raDec2TP, tanWcs);
TanPix2RaDec CTP2RaDec(identity, CommonTangentPoint);
CTP2TP = GtransfoCompose(&raDec2TP, &CTP2RaDec);

// jump from one TP to an other:
TanRaDec2Pix raDec2CTP(identity, CommonTangentPoint);
// TanPix2RaDec TP2RaDec(identity, tanWcs->TangentPoint());
// TP2CTP = GtransfoCompose(&raDec2CTP, &TP2RaDec);
TanPix2RaDec TP2RaDec(identity, tanWcs->TangentPoint());
TP2CTP = GtransfoCompose(&raDec2CTP, &TP2RaDec);
sky2TP = new TanRaDec2Pix(identity, tanWcs->TangentPoint());

// this one is needed for matches :
pix2CommonTangentPlane = GtransfoCompose(&raDec2CTP, tanWcs);

boresightRaDec = visitInfo->getBoresightRaDec();
airMass = visitInfo->getBoresightAirmass();
mjd = visitInfo->getDate().get(lsst::daf::base::DateTime::MJD);
double latitude = visitInfo->getObservatory().getLatitude();
double lst_obs = visitInfo->getEra();
double ra = visitInfo->getBoresightRaDec().getRa();
double dec = visitInfo->getBoresightRaDec().getDec();
double hourAngle = visitInfo->getBoresightHourAngle();
airMass = visitInfo->getBoresightAirmass();
mjd = visitInfo->getDate().get(lsst::daf::base::DateTime::MJD);

// lsstSim doesn't manage ERA (and thus Hour Angle) properly, so it's going to be NaN.
// Because we need the refraction vector later, go with 0 HA to prevent crashes on that NaN.
Expand All @@ -139,9 +110,36 @@ CcdImage::CcdImage(lsst::afw::table::SortedCatalogT<lsst::afw::table::SourceReco
tgz = sinz/cosz;
sineta = cos(latitude)*sin(hourAngle)/sinz;
coseta = sqrt(1 - sineta*sineta);
if (dec > latitude) coseta = -coseta;
if (boresightRaDec.getDec() > latitude) coseta = -coseta;
}
}

void CcdImage::setCommonTangentPoint(const Point &commonTangentPoint)
{
_commonTangentPoint = commonTangentPoint;

// use some other variable in case we later have to actually convert the
// as-read wcs:
const BaseTanWcs* tanWcs = readWcs.get();

/* we don't assume here that we know the internals of TanPix2RaDec:
to construct pix->TP, we do pix->sky->TP, although pix->sky
actually goes through TP */
GtransfoLin identity;
TanRaDec2Pix raDec2TP(identity, tanWcs->TangentPoint());
pix2TP = GtransfoCompose(&raDec2TP, tanWcs);
TanPix2RaDec CTP2RaDec(identity, commonTangentPoint);
CTP2TP = GtransfoCompose(&raDec2TP, &CTP2RaDec);

// jump from one TP to an other:
TanRaDec2Pix raDec2CTP(identity, commonTangentPoint);
TanPix2RaDec TP2RaDec(identity, tanWcs->TangentPoint());
TP2CTP = GtransfoCompose(&raDec2CTP, &TP2RaDec);
sky2TP = new TanRaDec2Pix(identity, tanWcs->TangentPoint());

// this one is needed for matches :
pix2CommonTangentPlane = GtransfoCompose(&raDec2CTP, tanWcs);
}

}
} // end of namespaces

0 comments on commit b0ecb53

Please sign in to comment.