Skip to content

Commit

Permalink
CRVAL is a Coord; pixel scale is an Angle; make SourceMatch use Angle…
Browse files Browse the repository at this point in the history 
… also (mostly)
  • Loading branch information
dstn committed May 10, 2011
1 parent f7c2b9d commit 2b49698
Show file tree
Hide file tree
Showing 16 changed files with 150 additions and 173 deletions.
6 changes: 4 additions & 2 deletions include/lsst/afw/detection/SourceMatch.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@
#include "lsst/afw/detection/Source.h"
#include "lsst/daf/base/Persistable.h"
#include "lsst/daf/base/PropertySet.h"
#include "lsst/afw/geom/Angle.h"

// forward declaration
namespace lsst { namespace afw { namespace formatters {
Expand All @@ -47,6 +48,7 @@ namespace lsst { namespace afw { namespace detection {
struct SourceMatch {
Source::Ptr first;
Source::Ptr second;
// match distance, in RADIANS or PIXELS depending on the type of match requested.
double distance;

SourceMatch() : first(), second(), distance(0.0) {}
Expand All @@ -56,8 +58,8 @@ struct SourceMatch {
};

std::vector<SourceMatch> matchRaDec(SourceSet const &set1, SourceSet const &set2,
double radius, bool closest=true);
std::vector<SourceMatch> matchRaDec(SourceSet const &set, double radius, bool symmetric = true);
lsst::afw::geom::Angle radius, bool closest=true);
std::vector<SourceMatch> matchRaDec(SourceSet const &set, lsst::afw::geom::Angle radius, bool symmetric = true);
std::vector<SourceMatch> matchXy(SourceSet const &set1, SourceSet const &set2,
double radius, bool closest=true);
std::vector<SourceMatch> matchXy(SourceSet const &set, double radius, bool symmetric = true);
Expand Down
10 changes: 9 additions & 1 deletion include/lsst/afw/geom/Angle.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,8 +21,9 @@ double const PI = boost::math::constants::pi<double>(); ///< The ratio of a circ
double const TWOPI = boost::math::constants::pi<double>() * 2.0;
double const HALFPI = boost::math::constants::pi<double>() * 0.5;

// NOTE, if you add things here, remember to also add them to
// NOTE, if you add things here, you must also add them to
// python/lsst/afw/geom/__init__.py
// (if you want them to accessible from python)

#if 0 && !defined(M_PI) // a good idea, but with ramifications
# define M_PI ::lsst::afw::geom::PI
Expand Down Expand Up @@ -96,6 +97,13 @@ class Angle {
/** Return an Angle's value as a double in arcseconds */
double asArcseconds() const { return asAngularUnits(arcseconds); }

double toUnitSphereDistanceSquared() const { return 2. * (1. - std::cos(asRadians())); }
// == 4.0 * pow(std::sin(0.5 * asRadians()), 2.0)
static Angle fromUnitSphereDistanceSquared(double d2) {
return (std::acos(1. - d2/2.)) * radians;
// == 2.0 * asin(0.5 * sqrt(d2))
}

/** Wraps this angle to the range [0, 2 pi) */
void wrap() {
_val = std::fmod(_val, TWOPI);
Expand Down
4 changes: 2 additions & 2 deletions include/lsst/afw/image/TanWcs.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,8 +90,8 @@ namespace image {

bool operator==(const TanWcs &) const;

// Returns the pixel scale, in arcsec/pixel.
double pixelScale() const;
// Returns the pixel scale, in Angle/pixel.
lsst::afw::geom::Angle pixelScale() const;

// Applies the SIP AP and BP distortion (used in the skyToPixel direction)
lsst::afw::geom::Point2D distortPixel(const lsst::afw::geom::Point2D pixel) const;
Expand Down
78 changes: 45 additions & 33 deletions include/lsst/afw/image/Wcs.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,15 +54,19 @@ namespace image {
/// - Greisen & Calabretta, 2002 A&A 395, 1061
/// - Calabretta & Greisen, 2002, A&A 395, 1077
///
/// In it's simplest sense, Wcs is used to convert from position in the sky (in right ascension
/// and declination) to pixel position on an image (and back again). It is, however, much more general
/// than that and can understand a myriad of different coordinate systems.
/// In its simplest sense, Wcs is used to convert from position in the sky (in
/// right ascension and declination) to pixel position on an image (and back
/// again). It is, however, much more general than that and can understand a
/// myriad of different coordinate systems.
///
/// A wcs can be constructed from a reference position (crval, crpix) and a translation matrix. Alternatively,
/// if you have the header from a fits file, you can create a Wcs object with the makeWcs() function. This
/// function determines whether your Wcs is one the subset of projection systems that is dealt with specially
/// by Lsst, and creates an object of the correct class. Otherwise, a pointer to a Wcs object is returned.
/// Most astronomical images use tangent plane projection, so makeWcs() returns a TanWcs object pointer
/// A wcs can be constructed from a reference position (crval, crpix) and a
/// translation matrix. Alternatively, if you have the header from a fits file,
/// you can create a Wcs object with the makeWcs() function. This function
/// determines whether your Wcs is one the subset of projection systems that is
/// dealt with specially by Lsst, and creates an object of the correct
/// class. Otherwise, a pointer to a Wcs object is returned. Most astronomical
/// images use tangent plane projection, so makeWcs() returns a TanWcs object
/// pointer
///
/// \code
/// import lsst.afw.image as afwImg
Expand All @@ -83,11 +87,12 @@ namespace image {
/// This class is implemented in by calls to the wcslib library
/// by Mark Calabretta http://www.atnf.csiro.au/people/mcalabre/WCS/
///
/// Note that we violate the Wcs standard in one minor way. The standard states that none
/// of the CRPIX or CRVAL keywords are required, for the header to be valid, and the appropriate values
/// should be set to 0.0 if the keywords are absent. This is a recipe for painful bugs in analysis, so
/// we violate the standard by insisting that the keywords CRPIX[1,2] and CRVAL[1,2] are present when
/// reading a header (keywords CRPIX1a etc are also accepted)
/// Note that we violate the Wcs standard in one minor way. The standard states
/// that none of the CRPIX or CRVAL keywords are required, for the header to be
/// valid, and the appropriate values should be set to 0.0 if the keywords are
/// absent. This is a recipe for painful bugs in analysis, so we violate the
/// standard by insisting that the keywords CRPIX[1,2] and CRVAL[1,2] are
/// present when reading a header (keywords CRPIX1a etc are also accepted)

class Wcs : public lsst::daf::base::Persistable,
public lsst::daf::data::LsstBase
Expand All @@ -96,7 +101,6 @@ class Wcs : public lsst::daf::base::Persistable,
typedef boost::shared_ptr<lsst::afw::image::Wcs> Ptr;
typedef boost::shared_ptr<lsst::afw::image::Wcs const> ConstPtr;

//Constructors
Wcs();
//Create a Wcs of the correct class using a fits header.
friend Wcs::Ptr makeWcs(lsst::daf::base::PropertySet::Ptr fitsMetadata,
Expand All @@ -113,38 +117,44 @@ class Wcs : public lsst::daf::base::Persistable,

bool operator==(const Wcs &) const;

//Accessors
lsst::afw::coord::Coord::Ptr getSkyOrigin() const; //Return crval
lsst::afw::geom::Point2D getPixelOrigin() const; //Return crpix
Eigen::Matrix2d getCDMatrix() const; //Return CD matrix
// Returns CRVAL
lsst::afw::coord::Coord::Ptr getSkyOrigin() const;
// Returns CRPIX
lsst::afw::geom::Point2D getPixelOrigin() const;
// Returns CD matrix. You would never have guessed that from the name.
Eigen::Matrix2d getCDMatrix() const;

virtual PTR(lsst::daf::base::PropertyList) getFitsMetadata() const;

/// Return true iff Wcs is valid
// Return true iff Wcs is valid
operator bool() const { return _nWcsInfo != 0; }

bool isFlipped() const; //Does the Wcs follow the convention of North=Up, East=Left or not

///Sky area covered by a pixel at position \c pix00 in units of square degrees.
// Does the Wcs follow the convention of North=Up, East=Left or not
// [This actually just measures the sign of the determinant of the CD matrix
// to determine the "handedness" of the coordinate system.]
bool isFlipped() const;

// Sky area covered by a pixel at position \c pix00 in units of square degrees.
double pixArea(lsst::afw::geom::Point2D pix00) const;

// Returns the pixel scale, in arcsec/pixel.
double pixelScale() const;
// Returns the pixel scale [Angle/pixel]
lsst::afw::geom::Angle pixelScale() const;

bool isInitialized() const;

//Convert from raDec to pixel space. Formerly called raDecToXY() and
//xyToRaDec(), but the name now reflects their increased generality. They may be
//used, e.g. to convert xy to Galactic coordinates
// Convert from raDec to pixel space. Formerly called raDecToXY() and
// xyToRaDec(), but the name now reflects their increased generality. They
// may be used, e.g. to convert xy to Galactic coordinates
lsst::afw::coord::Coord::Ptr pixelToSky(double pix1, double pix2) const;
lsst::afw::coord::Coord::Ptr pixelToSky(const lsst::afw::geom::Point2D pixel) const;

/// \note This routine is designed for the knowledgeable user in need of performance; it's safer to call
/// the version that returns a CoordPtr
/// \note This routine is designed for the knowledgeable user in need of
/// performance; it's safer to call the version that returns a CoordPtr
void pixelToSky(double pixel1, double pixel2, lsst::afw::geom::Angle& sky1, lsst::afw::geom::Angle& sky2) const;

// ASSUMES the angles are in the appropriate system for this WCS.
// ASSUMES the angles are in the appropriate coordinate system for this WCS.
lsst::afw::geom::Point2D skyToPixel(lsst::afw::geom::Angle sky1, lsst::afw::geom::Angle sky2) const;

lsst::afw::geom::Point2D skyToPixel(lsst::afw::coord::Coord::ConstPtr coord) const;
lsst::afw::geom::Point2D skyToIntermediateWorldCoord(lsst::afw::coord::Coord::ConstPtr coord) const;

Expand Down Expand Up @@ -227,7 +237,6 @@ class Wcs : public lsst::daf::base::Persistable,
int _wcshdrCtrl; ///< Controls messages to stderr from wcshdr (0 for none); see wcshdr.h for details
int _nReject;
lsst::afw::coord::CoordSystem _coordSystem;
bool _skyCoordsReversed;
};

namespace detail {
Expand All @@ -238,8 +247,11 @@ namespace detail {
}

Wcs::Ptr makeWcs(lsst::daf::base::PropertySet::Ptr fitsMetadata, bool stripMetadata=false);

Wcs::Ptr makeWcs(lsst::afw::geom::Point2D crval, lsst::afw::geom::Point2D crpix,

/*
Note, CD matrix elements must be in degrees/pixel.
*/
Wcs::Ptr makeWcs(lsst::afw::coord::Coord::ConstPtr crval, lsst::afw::geom::Point2D crpix,
double CD11, double CD12, double CD21, double CD22);

namespace detail {
Expand Down
19 changes: 0 additions & 19 deletions python/lsst/afw/image/imageLib.i
View file
Original file line number Diff line number Diff line change
Expand Up @@ -225,25 +225,6 @@ SWIG_SHARED_PTR_DERIVED(TanWcs, lsst::afw::image::Wcs, lsst::afw::image::TanWcs)
}
%}


%inline {
/**
* Create a WCS from crval, image, and the elements of CD
*/
lsst::afw::image::Wcs::Ptr createWcs(lsst::afw::geom::PointD crval,
lsst::afw::geom::PointD crpix,
double CD11, double CD12, double CD21, double CD22) {

Eigen::Matrix2d CD;
CD(0, 0) = CD11;
CD(0, 1) = CD12;
CD(1, 0) = CD21;
CD(1, 1) = CD22;

return lsst::afw::image::Wcs::Ptr(new lsst::afw::image::Wcs(crval, crpix, CD));
}
}

/************************************************************************************************************/

#if !defined(CAMERA_GEOM_LIB_I)
Expand Down

0 comments on commit 2b49698

Please sign in to comment.