/***************************************************************************
* GSkymap.cpp - Class that implements a sky map *
* ----------------------------------------------------------------------- *
* copyright (C) 2010-2012 by Juergen Knoedlseder *
* ----------------------------------------------------------------------- *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see . *
* *
***************************************************************************/
/**
* @file GSkymap.cpp
* @brief Sky map class implementation
* @author Juergen Knoedlseder
*/
/* __ Includes ___________________________________________________________ */
#ifdef HAVE_CONFIG_H
#include
#endif
#include "GException.hpp"
#include "GTools.hpp"
#include "GSkymap.hpp"
#include "GWcsRegistry.hpp"
#include "GWcslib.hpp"
#include "GWcsHPX.hpp"
#include "GFits.hpp"
#include "GFitsTableDoubleCol.hpp"
#include "GFitsImageDouble.hpp"
/* __ Method name definitions ____________________________________________ */
#define G_CONSTRUCT_HPX "GSkymap::GSkymap(std::string,std::string,int," \
"std::string,int)"
#define G_CONSTRUCT_MAP "GSkymap::GSkymap(std::string,std::string,GSkyDir," \
"int,int,double,double,int)"
#define G_OP_ACCESS_1D "GSkymap::operator(int&,int&)"
#define G_OP_ACCESS_2D "GSkymap::operator(GSkyPixel&,int&)"
#define G_OP_VALUE "GSkymap::operator(GSkyDir&,int&)"
#define G_READ "GSkymap::read(const GFitsHDU*)"
#define G_PIX2DIR "GSkymap::pix2dir(int)"
#define G_DIR2PIX "GSkymap::dir2pix(GSkyDir)"
#define G_XY2DIR "GSkymap::xy2dir(GSkyPixel)"
#define G_DIR2XY "GSkymap::dir2xy(GSkyDir)"
#define G_OMEGA1 "GSkymap::omega(int)"
#define G_OMEGA2 "GSkymap::omega(GSkyPixel)"
#define G_SET_WCS "GSkymap::set_wcs(std::string,std::string,double,double," \
"double,double,double,double,GMatrix,GVector)"
#define G_READ_HEALPIX "GSkymap::read_healpix(GFitsTable*)"
#define G_READ_WCS "GSkymap::read_wcs(GFitsImage*)"
#define G_ALLOC_WCS "GSkymap::alloc_wcs(GFitsImage*)"
/* __ Macros _____________________________________________________________ */
/* __ Coding definitions _________________________________________________ */
/* __ Debug definitions __________________________________________________ */
//#define G_READ_HEALPIX_DEBUG // Debug read_healpix
//#define G_READ_WCS_DEBUG // Debug read_wcs
/* __ Prototype __________________________________________________________ */
/*==========================================================================
= =
= Constructors/destructors =
= =
==========================================================================*/
/***********************************************************************//**
* @brief Void constructor
***************************************************************************/
GSkymap::GSkymap(void)
{
// Initialise class members for clean destruction
init_members();
// Return
return;
}
/***********************************************************************//**
* @brief FITS file constructor
*
* @param[in] filename FITS file name.
***************************************************************************/
GSkymap::GSkymap(const std::string& filename)
{
// Initialise class members for clean destruction
init_members();
// Load skymap
load(filename);
// Return
return;
}
/***********************************************************************//**
* @brief Healpix constructor
*
* @param[in] wcs World Coordinate System (HPX).
* @param[in] coords Coordinate System (CEL or GAL).
* @param[in] nside Nside parameter.
* @param[in] order Pixel ordering (RING or NEST).
* @param[in] nmaps Number of maps in set (default=1).
*
* @exception GException::wcs_invalid
* Invalid wcs parameter.
* @exception GException::skymap_bad_par
* Invalid sky map parameter.
***************************************************************************/
GSkymap::GSkymap(const std::string& wcs, const std::string& coords,
const int& nside, const std::string& order,
const int nmaps)
{
// Initialise class members for clean destruction
init_members();
// Check if wcs is HPX
if (toupper(wcs) != "HPX") {
throw GException::wcs_invalid(G_CONSTRUCT_HPX, wcs,
"WCS parameter must be 'HPX'.");
}
// Check if nmaps parameter is >0
if (nmaps < 1) {
throw GException::skymap_bad_par(G_CONSTRUCT_HPX, nmaps,
"nmaps parameter must be >0.");
}
// Allocate WCS
m_wcs = new GWcsHPX(nside, order, coords);
// Set number of pixels and number of maps
m_num_pixels = static_cast(m_wcs)->npix();
m_num_maps = nmaps;
// Allocate pixels
alloc_pixels();
// Return
return;
}
/***********************************************************************//**
* @brief WCS sky map constructor
*
* @param[in] wcs World Coordinate System.
* @param[in] coords Coordinate System (CEL or GAL).
* @param[in] x X coordinate of sky map centre (deg).
* @param[in] y Y coordinate of sky map centre (deg).
* @param[in] dx Pixel size in x direction at centre (deg/pixel).
* @param[in] dy Pixel size in y direction at centre (deg/pixel).
* @param[in] nx Number of pixels in x direction.
* @param[in] ny Number of pixels in y direction.
* @param[in] nmaps Number of maps in set (default=1).
*
* @exception GException::skymap_bad_par
* Invalid sky map parameter.
***************************************************************************/
GSkymap::GSkymap(const std::string& wcs, const std::string& coords,
double const& x, double const& y,
double const& dx, double const& dy,
const int& nx, const int& ny, const int nmaps)
{
// Initialise class members for clean destruction
init_members();
// Check parameters
if (nx < 1) {
throw GException::skymap_bad_par(G_CONSTRUCT_MAP, nx,
"nx parameter must be >0.");
}
if (ny < 1) {
throw GException::skymap_bad_par(G_CONSTRUCT_MAP, ny,
"ny parameter must be >0.");
}
if (nmaps < 1) {
throw GException::skymap_bad_par(G_CONSTRUCT_MAP, nmaps,
"nmaps parameter must be >0.");
}
// Set WCS
double crval1 = x;
double crval2 = y;
double crpix1 = double(nx+1)/2.0;
double crpix2 = double(ny+1)/2.0;
double cdelt1 = dx;
double cdelt2 = dy;
GMatrix cd(2,2);
GVector pv2(21);
set_wcs(wcs, coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2,
cd, pv2);
// Set number of pixels and number of maps
m_num_x = nx;
m_num_y = ny;
m_num_pixels = m_num_x * m_num_y;
m_num_maps = nmaps;
// Allocate pixels
alloc_pixels();
// Return
return;
}
/***********************************************************************//**
* @brief Copy constructor
*
* @param[in] map Sky map.
***************************************************************************/
GSkymap::GSkymap(const GSkymap& map)
{
// Initialise class members for clean destruction
init_members();
// Copy members
copy_members(map);
// Return
return;
}
/***********************************************************************//**
* @brief Destructor
***************************************************************************/
GSkymap::~GSkymap(void)
{
// Free members
free_members();
// Return
return;
}
/*==========================================================================
= =
= Operators =
= =
==========================================================================*/
/***********************************************************************//**
* @brief 1D pixel access operator
*
* @param[in] pixel Pixel index (0,1,...,m_num_pixels).
* @param[in] map Map index (0,1,...,m_num_maps).
*
* @exception GException::out_of_range
* Pixel index and/or map index are outside valid range.
*
* Access sky map pixel by its index, where the most quickly varying axis is
* the x axis of the map.
*
* @todo Implement proper skymap exception (actual is for matrix elements)
***************************************************************************/
double& GSkymap::operator() (const int& pixel, const int& map)
{
// Throw an error if pixel index or map index is not in valid range
#if defined(G_RANGE_CHECK)
if (pixel < 0 || pixel >= m_num_pixels) {
throw GException::out_of_range(G_OP_ACCESS_1D, pixel, 0, m_num_pixels-1);
}
if (map < 0 || map >= m_num_maps) {
throw GException::out_of_range(G_OP_ACCESS_1D, map, 0, m_num_maps-1);
}
#endif
// Return reference to pixel value
return m_pixels[pixel+m_num_pixels*map];
}
/***********************************************************************//**
* @brief 1D pixel access operator (const variant)
*
* @param[in] pixel Pixel index (0,1,...,m_num_pixels).
* @param[in] map Map index (0,1,...,m_num_maps).
*
* @exception GException::out_of_range
* Pixel index and/or map index are outside valid range.
*
* Access sky map pixel by its index, where the most quickly varying axis is
* the x axis of the map.
*
* @todo Implement proper skymap exception (actual is for matrix elements)
***************************************************************************/
const double& GSkymap::operator() (const int& pixel, const int& map) const
{
// Throw an error if pixel index or map index is not in valid range
#if defined(G_RANGE_CHECK)
if (pixel < 0 || pixel >= m_num_pixels) {
throw GException::out_of_range(G_OP_ACCESS_1D, pixel, 0, m_num_pixels-1);
}
if (map < 0 || map >= m_num_maps) {
throw GException::out_of_range(G_OP_ACCESS_1D, map, 0, m_num_maps-1);
}
#endif
// Return reference to pixel value
return m_pixels[pixel+m_num_pixels*map];
}
/***********************************************************************//**
* @brief 2D pixel access operator
*
* @param[in] pixel Sky pixel.
* @param[in] map Map index (0,1,...,m_num_maps).
*
* @exception GException::out_of_range
* Sky pixel and/or map index are outside valid range.
*
* Access sky map pixel by its 2D index (x,y) that is implemented by the
* GSkyPixel class.
*
* @todo Implement proper skymap exception (actual is for matrix elements)
***************************************************************************/
double& GSkymap::operator() (const GSkyPixel& pixel, const int& map)
{
// Throw an error if pixel index or map index is not in valid range
#if defined(G_RANGE_CHECK)
if (!isinmap(pixel)) {
throw GException::out_of_range(G_OP_ACCESS_2D,
int(pixel.x()), int(pixel.y()),
m_num_x-1, m_num_y-1);
}
if (map < 0 || map >= m_num_maps) {
throw GException::out_of_range(G_OP_ACCESS_2D, map, 0, m_num_maps-1);
}
#endif
// Get pixel index
int index = xy2pix(pixel);
// Return reference to pixel value
return m_pixels[index+m_num_pixels*map];
}
/***********************************************************************//**
* @brief 2D pixel access operator
*
* @param[in] pixel Sky pixel.
* @param[in] map Map index (0,1,...,m_num_maps).
*
* @exception GException::out_of_range
* Sky pixel and/or map index are outside valid range.
*
* Access sky map pixel by its 2D index (x,y) that is implemented by the
* GSkyPixel class.
*
* @todo Implement proper skymap exception (actual is for matrix elements)
***************************************************************************/
const double& GSkymap::operator() (const GSkyPixel& pixel, const int& map) const
{
// Throw an error if pixel index or map index is not in valid range
#if defined(G_RANGE_CHECK)
if (!isinmap(pixel)) {
throw GException::out_of_range(G_OP_ACCESS_2D,
int(pixel.x()), int(pixel.y()),
m_num_x-1, m_num_y-1);
}
if (map < 0 || map >= m_num_maps) {
throw GException::out_of_range(G_OP_ACCESS_2D, map, 0, m_num_maps-1);
}
#endif
// Get pixel index
int index = xy2pix(pixel);
// Return reference to pixel value
return m_pixels[index+m_num_pixels*map];
}
/***********************************************************************//**
* @brief Return interpolated skymap value for sky direction
*
* @param[in] dir Sky direction.
* @param[in] map Map index (0,1,...,m_num_maps).
*
* @exception GException::out_of_range
* Map index lies outside valid range.
*
* Returns the skymap value for a given sky direction, obtained by bi-linear
* interpolation of the neighbouring pixels. If the sky direction falls
* outside the area covered by the skymap, a value of 0 is returned.
***************************************************************************/
double GSkymap::operator() (const GSkyDir& dir, const int& map) const
{
// Throw an error if the map index is not in valid range
#if defined(G_RANGE_CHECK)
if (map < 0 || map >= m_num_maps) {
throw GException::out_of_range(G_OP_VALUE, map, 0, m_num_maps-1);
}
#endif
// Initialise intensity
double intensity = 0.0;
// Determine sky pixel
GSkyPixel pixel = dir2xy(dir);
// Continue only if pixel is within the map
if (isinmap(pixel)) {
// Set left indices for interpolation. The left index is comprised
// between 0 and npixels-2. By definition, the right index is then
// the left index + 1
int inx_x = int(pixel.x());
int inx_y = int(pixel.y());
if (inx_x < 0) {
inx_x = 0;
}
else if (inx_x > m_num_x-2) {
inx_x = m_num_x - 2;
}
if (inx_y < 0) {
inx_y = 0;
}
else if (inx_y > m_num_y-2) {
inx_y = m_num_y - 2;
}
// Set weighting factors for interpolation
double wgt_x_right = (pixel.x() - inx_x);
double wgt_x_left = 1.0 - wgt_x_right;
double wgt_y_right = (pixel.y() - inx_y);
double wgt_y_left = 1.0 - wgt_y_right;
// Compute skymap pixel indices for bi-linear interpolation
int inx1 = inx_x + inx_y * m_num_x;
int inx2 = inx1 + m_num_x;
int inx3 = inx1 + 1;
int inx4 = inx2 + 1;
// Compute weighting factors for bi-linear interpolation
double wgt1 = wgt_x_left * wgt_y_left;
double wgt2 = wgt_x_left * wgt_y_right;
double wgt3 = wgt_x_right * wgt_y_left;
double wgt4 = wgt_x_right * wgt_y_right;
// Compute map offset
int offset = m_num_pixels * map;
// Compute interpolated skymap value
intensity = wgt1 * m_pixels[inx1 + offset] +
wgt2 * m_pixels[inx2 + offset] +
wgt3 * m_pixels[inx3 + offset] +
wgt4 * m_pixels[inx4 + offset];
} // endif: pixel was within map
// Return intensity
return intensity;
}
/***********************************************************************//**
* @brief Assignment operator
*
* @param[in] map Sky map.
***************************************************************************/
GSkymap& GSkymap::operator= (const GSkymap& map)
{
// Execute only if object is not identical
if (this != &map) {
// Free members
free_members();
// Initialise private members for clean destruction
init_members();
// Copy members
copy_members(map);
} // endif: object was not identical
// Return this object
return *this;
}
/*==========================================================================
= =
= Public methods =
= =
==========================================================================*/
/***********************************************************************//**
* @brief Clear instance.
*
* This method properly resets the object to an initial state.
***************************************************************************/
void GSkymap::clear(void)
{
// Free class members
free_members();
// Initialise members
init_members();
// Return
return;
}
/***********************************************************************//**
* @brief Clone instance
***************************************************************************/
GSkymap* GSkymap::clone(void) const
{
return new GSkymap(*this);
}
/***********************************************************************//**
* @brief Load skymap from FITS file.
*
* @param[in] filename FITS file name..
*
* Loads HEALPix and non HEALPix skymaps. First searches for HEALPix map in
* FITS file by scanning all HDUs for PIXTYPE=HEALPIX. If no HEALPix map has
* been found then search load first non-empty image.
***************************************************************************/
void GSkymap::load(const std::string& filename)
{
// Free memory and initialise members
free_members();
init_members();
// Open FITS file
GFits fits(filename);
// Get number of HDUs
int num = fits.size();
// Initialize pointer to HDU and load flag
GFitsHDU* hdu = NULL;
bool loaded = false;
// First search for HEALPix extension. We can skip the first extension
// since this is always an image and a HEALPix map is stored in a
// binary table
for (int extno = 1; extno < num; ++extno) {
// Get pointer to HDU
hdu = fits.hdu(extno);
// If PIXTYPE keyword equals "HEALPIX" then load map
try {
if (hdu->string("PIXTYPE") == "HEALPIX") {
read_healpix(static_cast(hdu));
loaded = true;
break;
}
}
catch (GException::fits_key_not_found &e) {
}
} // endfor: looped over HDUs
// If we have not found a HEALPIX map then search now for image.
// Skip empty images
if (!loaded) {
for (int extno = 0; extno < num; ++extno) {
// Get pointer to HDU
hdu = fits.hdu(extno);
// Skip if extension is not an image
if (extno > 0) {
if (hdu->string("XTENSION") != "IMAGE")
continue;
}
// Load WCS map
read_wcs(static_cast(hdu));
loaded = true;
break;
} // endfor: looped over HDUs
} // endif: no HEALPix map found
// Close FITS file
fits.close();
// Return
return;
}
/***********************************************************************//**
* @brief Save skymap into FITS file.
*
* @param[in] filename FITS file name.
* @param[in] clobber Overwrite existing file? (true=yes)
*
* The method does nothing if the skymap holds no valid WCS.
***************************************************************************/
void GSkymap::save(const std::string& filename, bool clobber) const
{
// Continue only if we have data to save
if (m_wcs != NULL) {
// Initialise HDU pointer
GFitsHDU* hdu = NULL;
// Case A: Skymap is Healpix
if (m_wcs->code() == "HPX") {
hdu = create_healpix_hdu();
}
// Case B: Skymap is not Healpix
else {
hdu = create_wcs_hdu();
}
// Create FITS file and save it to disk
if (hdu != NULL) {
GFits fits;
fits.append(*hdu);
fits.saveto(filename, clobber);
}
// Delete HDU
if (hdu != NULL) delete hdu;
} // endif: we had data to save
// Return
return;
}
/***********************************************************************//**
* @brief Read skymap from FITS HDU
*
* @param[in] hdu FITS HDU.
*
* The method returns an empty skymap of the HDU pointer was not valid.
***************************************************************************/
void GSkymap::read(const GFitsHDU* hdu)
{
// Free memory and initialise members
free_members();
init_members();
// Continue only if HDU pointer is valid
if (hdu != NULL) {
// Initialize load flag
bool loaded = false;
// Try load as HEALPix map
try {
if (hdu->string("PIXTYPE") == "HEALPIX") {
read_healpix(static_cast(hdu));
loaded = true;
}
}
catch (GException::fits_key_not_found &e) {
}
// ... otherwise try loading as non HEALPix map
if (!loaded) {
// Load only if HDU contains an image
if (hdu->exttype() == 0) {
read_wcs(static_cast(hdu));
loaded = true;
}
} // endif
} // endif: HDU pointer was valid
// Return
return;
}
/***********************************************************************//**
* @brief Write skymap into FITS file
*
* @param[in] file FITS file pointer.
***************************************************************************/
void GSkymap::write(GFits* file) const
{
// Continue only if we have data to save
if (m_wcs != NULL) {
// Initialise HDU pointer
GFitsHDU* hdu = NULL;
// Case A: Skymap is Healpix
if (m_wcs->code() == "HPX") {
hdu = create_healpix_hdu();
}
// Case B: Skymap is not Healpix
else {
hdu = create_wcs_hdu();
}
// Append HDU to FITS file.
if (hdu != NULL) {
file->append(*hdu);
}
// Delete HDU
if (hdu != NULL) delete hdu;
} // endif: we had data to save
// Return
return;
}
/***********************************************************************//**
* @brief Returns sky direction of pixel
*
* @param[in] pix Pixel number (0,1,...,m_num_pixels).
*
* @exception GException::wcs
* No valid WCS found.
*
* Returns sky direction for a given sky map pixel. This methods works for
* sky maps that have a 1D (e.g. HEALPix) or 2D pixel indexation scheme.
* Automatic index conversion is provided so that 1D schemes need only to
* implement 1D conversion methods while 2D schemes need only to implement
* 2D conversion methods.
***************************************************************************/
GSkyDir GSkymap::pix2dir(const int& pix) const
{
// Throw error if WCS is not valid
if (m_wcs == NULL) {
throw GException::wcs(G_PIX2DIR, "No valid WCS found.");
}
// Determine sky direction from pixel. Use 2D version if sky map is
// 2D, otherwise use 1D version.
GSkyDir dir = (m_num_x == 0) ? m_wcs->pix2dir(pix)
: m_wcs->xy2dir(pix2xy(pix));
// Return sky direction
return dir;
}
/***********************************************************************//**
* @brief Returns pixel index for a given sky direction
*
* @param[in] dir Sky direction.
*
* @exception GException::wcs
* No valid WCS found.
*
* Returns sky map pixel for a given sky direction. This methods works for
* sky maps that have a 1D (e.g. HEALPix) or 2D pixel indexation scheme.
* Automatic index conversion is provided so that 1D schemes need only to
* implement 1D conversion methods while 2D schemes need only to implement
* 2D conversion methods.
***************************************************************************/
int GSkymap::dir2pix(const GSkyDir& dir) const
{
// Throw error if WCS is not valid
if (m_wcs == NULL) {
throw GException::wcs(G_DIR2PIX, "No valid WCS found.");
}
// Determine 1D pixel index for a given sky direction
int pix = (m_num_x == 0) ? m_wcs->dir2pix(dir)
: xy2pix(m_wcs->dir2xy(dir));
// Return pixel index
return pix;
}
/***********************************************************************//**
* @brief Returns sky direction of pixel
*
* @param[in] pix Sky map pixel.
*
* @exception GException::wcs
* No valid WCS found.
*
* Returns sky direction for a given sky map pixel. This methods works for
* sky maps that have a 1D (e.g. HEALPix) or 2D pixel indexation scheme.
* Automatic index conversion is provided so that 1D schemes need only to
* implement 1D conversion methods while 2D schemes need only to implement
* 2D conversion methods.
***************************************************************************/
GSkyDir GSkymap::xy2dir(const GSkyPixel& pix) const
{
// Throw error if WCS is not valid
if (m_wcs == NULL) {
throw GException::wcs(G_XY2DIR, "No valid WCS found.");
}
// Determine sky direction from pixel. Use 2D version if sky map is
// 2D, otherwise use 1D version.
GSkyDir dir = (m_num_x == 0) ? m_wcs->pix2dir(xy2pix(pix))
: m_wcs->xy2dir(pix);
// Return sky direction
return dir;
}
/***********************************************************************//**
* @brief Returns sky map pixel for a given sky direction
*
* @param[in] dir Sky direction.
*
* @exception GException::wcs
* No valid WCS found.
*
* Returns sky map pixel for a given sky direction. This methods works for
* sky maps that have a 1D (e.g. HEALPix) or 2D pixel indexation scheme.
* Automatic index conversion is provided so that 1D schemes need only to
* implement 1D conversion methods while 2D schemes need only to implement
* 2D conversion methods.
***************************************************************************/
GSkyPixel GSkymap::dir2xy(const GSkyDir& dir) const
{
// Throw error if WCS is not valid
if (m_wcs == NULL) {
throw GException::wcs(G_DIR2XY, "No valid WCS found.");
}
// Determine 1D pixel index for a given sky direction
GSkyPixel pix = (m_num_x == 0) ? pix2xy(m_wcs->dir2pix(dir))
: m_wcs->dir2xy(dir);
// Return pixel index
return pix;
}
/***********************************************************************//**
* @brief Returns solid angle of pixel
*
* @param[in] pix Pixel number (0,1,...,m_num_pixels).
*
* @exception GException::wcs
* No valid WCS found.
***************************************************************************/
double GSkymap::omega(const int& pix) const
{
// Throw error if WCS is not valid
if (m_wcs == NULL) {
throw GException::wcs(G_OMEGA1, "No valid WCS found.");
}
// Determine solid angle from pixel. Use 2D version if sky map is
// 2D, otherwise use 1D version.
double omega = (m_num_x == 0) ? m_wcs->omega(pix)
: m_wcs->omega(pix2xy(pix));
// Return solid angle
return omega;
}
/***********************************************************************//**
* @brief Returns solid angle of pixel
*
* @param[in] pix Sky map pixel.
*
* @exception GException::wcs
* No valid WCS found.
***************************************************************************/
double GSkymap::omega(const GSkyPixel& pix) const
{
// Throw error if WCS is not valid
if (m_wcs == NULL) {
throw GException::wcs(G_OMEGA2, "No valid WCS found.");
}
// Determine solid angle from pixel. Use 2D version if sky map is
// 2D, otherwise use 1D version.
double omega = (m_num_x == 0) ? m_wcs->omega(xy2pix(pix))
: m_wcs->omega(pix);
// Return solid angle
return omega;
}
/***********************************************************************//**
* @brief Returns number of pixels
***************************************************************************/
int GSkymap::npix(void) const
{
// Return number of pixels
return m_num_pixels;
}
/***********************************************************************//**
* @brief Returns number of pixels in x coordinate
***************************************************************************/
int GSkymap::nx(void) const
{
// Return number of pixels
return m_num_x;
}
/***********************************************************************//**
* @brief Returns number of pixels in y coordinate
***************************************************************************/
int GSkymap::ny(void) const
{
// Return number of pixels
return m_num_y;
}
/***********************************************************************//**
* @brief Returns number of maps
***************************************************************************/
int GSkymap::nmaps(void) const
{
// Return number of maps
return m_num_maps;
}
/***********************************************************************//**
* @brief Converts 2D index (x,y) into 1D pixel index
*
* @param[in] pix 2D pixel index.
*
* The (x,y) value is rounded to nearest integers before conversion. The x
* axis is assumed as the most rapidely varying index.
***************************************************************************/
int GSkymap::xy2pix(const GSkyPixel& pix) const
{
// Get x and y indices by rounding the (x,y) values
int ix = int(pix.x()+0.5);
int iy = int(pix.y()+0.5);
// Return index
return (ix+iy*m_num_x);
}
/***********************************************************************//**
* @brief Converts 1D pixel index into 2D index (x,y)
*
* @param[in] pix 1D pixel index.
*
* If skymap is not 2D the pixel index is returned in the x element of the
* GSkyPixel object.
***************************************************************************/
GSkyPixel GSkymap::pix2xy(const int& pix) const
{
// Get x and y indices
double x;
double y;
if (m_num_x != 0) {
x = double(pix % m_num_x);
y = double(pix / m_num_x);
}
else {
x = double(pix);
y = 0.0;
}
// Set pixel
GSkyPixel pixel(x, y);
// Return pixel
return pixel;
}
/***********************************************************************//**
* @brief Set WCS skymap
*
* @param[in] wcs World Coordinate System.
*
* Sets the World Coordinate System of sky map. The method performs a deep
* copy of the input argument, allowing to destroy the argument after using
* the method.
***************************************************************************/
void GSkymap::wcs(const GWcs& wcs)
{
// Free any existing WCS
if (m_wcs != NULL) delete m_wcs;
// Clone input WCS
m_wcs = wcs.clone();
// Return
return;
}
/***********************************************************************//**
* @brief Verifies if sky direction falls in map
*
* @param[in] dir Sky direction.
*
* This method checks if the specified sky direction falls within the pixels
* covered by the skymap. The method uses the dir2xy method to convert the
* sky direction into 2D pixel indices, and then checks whether the pixel
* indices fall in the skymap.
***************************************************************************/
bool GSkymap::isinmap(const GSkyDir& dir) const
{
// Convert sky direction into sky pixel
GSkyPixel pixel = dir2xy(dir);
// Return location flag
return (isinmap(pixel));
}
/***********************************************************************//**
* @brief Verifies if sky pixel falls in map
*
* @param[in] pixel Sky pixel.
*
* This method checks if the specified sky pixel is within the skymap.
***************************************************************************/
bool GSkymap::isinmap(const GSkyPixel& pixel) const
{
// Initialise location flag
bool inmap = false;
// If pixel is in range then set location flag to true
if ((pixel.x()+0.5 >= 0.0 && pixel.x()-0.5 < m_num_x) &&
(pixel.y()+0.5 >= 0.0 && pixel.y()-0.5 < m_num_y)) {
inmap = true;
}
// Return location flag
return inmap;
}
/***********************************************************************//**
* @brief Print models
***************************************************************************/
std::string GSkymap::print(void) const
{
// Initialise result string
std::string result;
// Append header
result.append("=== GSkymap ===\n");
result.append(parformat("Number of pixels")+str(m_num_pixels)+"\n");
result.append(parformat("Number of maps")+str(m_num_maps));
if (m_wcs != NULL && m_wcs->code() != "HPX") {
result.append("\n"+parformat("X axis dimension")+str(m_num_x));
result.append("\n"+parformat("Y axis dimension")+str(m_num_y));
}
// Append WCS information
if (m_wcs != NULL) {
result.append("\n"+m_wcs->print());
}
else {
result.append("\n"+parformat("WCS")+"not defined");
}
// Return result
return result;
}
/*==========================================================================
= =
= Private methods =
= =
==========================================================================*/
/***********************************************************************//**
* @brief Initialise class members
***************************************************************************/
void GSkymap::init_members(void)
{
// Initialise members
m_num_pixels = 0;
m_num_maps = 0;
m_num_x = 0;
m_num_y = 0;
m_wcs = NULL;
m_pixels = NULL;
// Return
return;
}
/***********************************************************************//**
* @brief Allocate skymap pixels
***************************************************************************/
void GSkymap::alloc_pixels(void)
{
// Compute data size
int size = m_num_pixels * m_num_maps;
// Continue only if there are pixels
if (size > 0) {
// Allocate pixels and initialize them to 0
m_pixels = new double[size];
for (int i = 0; i < size; ++i) {
m_pixels[i] = 0.0;
}
} // endif: there were pixels
// Return
return;
}
/***********************************************************************//**
* @brief Copy class members
*
* @param[in] map Sky map.
***************************************************************************/
void GSkymap::copy_members(const GSkymap& map)
{
// Copy attributes
m_num_pixels = map.m_num_pixels;
m_num_maps = map.m_num_maps;
m_num_x = map.m_num_x;
m_num_y = map.m_num_y;
// Clone WCS if it is valid
if (map.m_wcs != NULL) m_wcs = map.m_wcs->clone();
// Compute data size
int size = m_num_pixels * m_num_maps;
// Copy pixels
if (size > 0 && map.m_pixels != NULL) {
alloc_pixels();
for (int i = 0; i < size; ++i) {
m_pixels[i] = map.m_pixels[i];
}
}
// Return
return;
}
/***********************************************************************//**
* @brief Delete class members
***************************************************************************/
void GSkymap::free_members(void)
{
// Free memory
if (m_wcs != NULL) delete m_wcs;
if (m_pixels != NULL) delete [] m_pixels;
// Signal free pointers
m_wcs = NULL;
m_pixels = NULL;
// Reset number of pixels
m_num_pixels = 0;
m_num_maps = 0;
m_num_x = 0;
m_num_y = 0;
// Return
return;
}
/***********************************************************************//**
* @brief Set WCS
*
* @param[in] wcs World Coordinate System code.
* @param[in] coords Coordinate system.
* @param[in] crval1 X value of reference pixel.
* @param[in] crval2 Y value of reference pixel.
* @param[in] crpix1 X index of reference pixel.
* @param[in] crpix2 Y index of reference pixel.
* @param[in] cdelt1 Increment in x direction at reference pixel (deg).
* @param[in] cdelt2 Increment in y direction at reference pixel (deg).
* @param[in] cd Astrometry parameters (2x2 matrix, deg/pixel).
* @param[in] pv2 Projection parameters (length WCS type dependent).
*
* @exception GException::wcs_invalid
* Invalid wcs parameter (World Coordinate System not known).
*
* This method sets the WCS projection pointer based on the WCS code and
* sky map parameters. It makes use of the GWcsRegistry class to allocate
* the correct derived class. Note that this method does not support the
* HPX projection.
*
* @todo Remove cd and pv2 parameters.
***************************************************************************/
void GSkymap::set_wcs(const std::string& wcs, const std::string& coords,
const double& crval1, const double& crval2,
const double& crpix1, const double& crpix2,
const double& cdelt1, const double& cdelt2,
const GMatrix& cd, const GVector& pv2)
{
// Convert WCS to upper case
std::string uwcs = toupper(wcs);
// Check if HPX was requested (since this is not allowed)
if (uwcs == "HPX") {
throw GException::wcs_invalid(G_SET_WCS, uwcs,
"Method not valid for HPX projection.");
}
// ... otherwise get projection from registry
else {
// Allocate WCS registry
GWcsRegistry registry;
// Allocate projection from registry
m_wcs = registry.alloc(uwcs);
// Signal if projection type is not known
if (m_wcs == NULL) {
std::string message = "Projection code not known. "
"Should be one of "+registry.list()+".";
throw GException::wcs_invalid(G_SET_WCS, uwcs, message);
}
// Setup WCS
static_cast(m_wcs)->set(coords, crval1, crval2, crpix1, crpix2,
cdelt1, cdelt2);
} // endelse: got projection from registry
// Return
return;
}
/***********************************************************************//**
* @brief Read Healpix data from FITS table.
*
* @param[in] hdu FITS HDU containing the Healpix data.
*
* HEALPix data may be stored in various formats depending on the
* application that has writted the data. HEALPix IDL, for example, may
* store the data in vectors of length 1024 if the number of pixels is
* a multiple of 1024. On the other hand, vectors may also be used to store
* several HEALPix maps into a single column. Alternatively, multiple maps
* may be stored in multiple columns.
***************************************************************************/
void GSkymap::read_healpix(const GFitsTable* hdu)
{
// Continue only if HDU is valid
if (hdu != NULL) {
// Determine number of rows and columns in table
int nrows = hdu->nrows();
int ncols = hdu->ncols();
#if defined(G_READ_HEALPIX_DEBUG)
std::cout << "nrows=" << nrows << " ncols=" << ncols << std::endl;
#endif
// Allocate Healpix WCS
m_wcs = new GWcsHPX;
// Read WCS information from FITS header
m_wcs->read(hdu);
// Set number of pixels based on NSIDE parameter
m_num_pixels = static_cast(m_wcs)->npix();
#if defined(G_READ_HEALPIX_DEBUG)
std::cout << "m_num_pixels=" << m_num_pixels << std::endl;
#endif
// Number of map pixels has to be a multiple of the number of
// rows in column
if (m_num_pixels % nrows != 0) {
throw GException::skymap_bad_size(G_READ_HEALPIX, nrows,
m_num_pixels);
}
// Determine vector length for HEALPix data storage
int nentry = m_num_pixels / nrows;
#if defined(G_READ_HEALPIX_DEBUG)
std::cout << "nentry=" << nentry << std::endl;
#endif
// Determine number of maps from the number of maps that fit into
// all columns. Only count columns that can fully hold the map.
m_num_maps = 0;
for (int icol = 0; icol < ncols; ++icol) {
const GFitsTableCol* col = &(*hdu)[icol];
if (col->number() % nentry == 0) {
m_num_maps += col->number() / nentry;
}
}
#if defined(G_READ_HEALPIX_DEBUG)
std::cout << "m_num_maps=" << m_num_maps << std::endl;
#endif
// Allocate pixels to hold the map
alloc_pixels();
// Initialise map counter
int imap = 0;
// Loop over all columns
for (int icol = 0; icol < ncols; ++icol) {
// Get next column
const GFitsTableCol* col = &(*hdu)[icol];
// Only consider columns that can fully hold maps
if (col->number() % nentry == 0) {
// Determine number of maps in column
int num = col->number() / nentry;
// Loop over all maps in column
int inx_start = 0;
int inx_end = nentry;
for (int i = 0; i < num; ++i) {
// Load map
double *ptr = m_pixels + m_num_pixels*imap;
for (int row = 0; row < col->length(); ++row) {
for (int inx = inx_start; inx < inx_end; ++inx) {
*ptr++ = col->real(row,inx);
}
}
#if defined(G_READ_HEALPIX_DEBUG)
std::cout << "Load map=" << imap << " index="
<< inx_start << "-" << inx_end << std::endl;
#endif
// Increment index range
inx_start = inx_end;
inx_end += nentry;
// Increment map counter
imap++;
// Break if we have loaded all maps
if (imap >= m_num_maps) {
break;
}
} // endfor: looped over all maps in column
} // endif: column could fully hold maps
// Break if we have loaded all maps
if (imap >= m_num_maps) {
break;
}
} // endfor: looped over all columns
} // endif: HDU was valid
// Return
return;
}
/***********************************************************************//**
* @brief Read WCS image from FITS HDU
*
* @param[in] hdu FITS HDU containing the WCS image.
*
* @exception GException::skymap_bad_image_dim
* WCS image has invalid dimension (naxis=2 or 3).
***************************************************************************/
void GSkymap::read_wcs(const GFitsImage* hdu)
{
// Continue only if HDU is valid
if (hdu != NULL) {
// Allocate WCS
alloc_wcs(hdu);
// Read WCS information from FITS header
m_wcs->read(hdu);
// Extract map dimension and number of maps from image
if (hdu->naxis() == 2) {
m_num_x = hdu->naxes(0);
m_num_y = hdu->naxes(1);
m_num_maps = 1;
}
else if (hdu->naxis() >= 3) {
m_num_x = hdu->naxes(0);
m_num_y = hdu->naxes(1);
m_num_maps = hdu->naxes(2);
}
else {
throw GException::skymap_bad_image_dim(G_READ_WCS, hdu->naxis());
}
#if defined(G_READ_WCS_DEBUG)
std::cout << "m_num_x=" << m_num_x << std::endl;
std::cout << "m_num_y=" << m_num_y << std::endl;
std::cout << "m_num_maps=" << m_num_maps << std::endl;
#endif
// Compute number of pixels
m_num_pixels = m_num_x * m_num_y;
#if defined(G_READ_WCS_DEBUG)
std::cout << "m_num_pixels=" << m_num_pixels << std::endl;
#endif
// Allocate pixels to hold the map
alloc_pixels();
// Read image
if (hdu->naxis() == 2) {
double* ptr = m_pixels;
for (int iy = 0; iy < m_num_y; ++iy) {
for (int ix = 0; ix < m_num_x; ++ix) {
*ptr++ = hdu->pixel(ix,iy);
}
}
}
else {
double* ptr = m_pixels;
for (int imap = 0; imap < m_num_maps; ++imap) {
for (int iy = 0; iy < m_num_y; ++iy) {
for (int ix = 0; ix < m_num_x; ++ix) {
*ptr++ = hdu->pixel(ix,iy,imap);
}
}
}
}
} // endif: HDU was valid
// Return
return;
}
/***********************************************************************//**
* @brief Allocate WCS class
*
* @param[in] hdu FITS HDU containing the WCS image.
*
* @exception GException::fits_key_not_found
* Unable to find required FITS header keyword.
* @exception GException::skymap_bad_ctype
* CTYPE1 and CTYPE2 keywords are incompatible.
* @exception GException::wcs_invalid
* WCS projection of FITS file not supported by GammaLib.
***************************************************************************/
void GSkymap::alloc_wcs(const GFitsImage* hdu)
{
// Continue only if HDU is valid
if (hdu != NULL) {
// Get standard keywords
std::string ctype1 = hdu->string("CTYPE1");
std::string ctype2 = hdu->string("CTYPE2");
// Extract projection type
std::string xproj = ctype1.substr(5,3);
std::string yproj = ctype2.substr(5,3);
// Check that projection type is identical on both axes
if (xproj != yproj) {
throw GException::skymap_bad_ctype(G_ALLOC_WCS,
ctype1, ctype2);
}
// Allocate WCS registry
GWcsRegistry registry;
// Allocate projection from registry
m_wcs = registry.alloc(xproj);
// Signal if projection type is not known
if (m_wcs == NULL) {
std::string message = "Projection code not known. "
"Should be one of "+registry.list()+".";
throw GException::wcs_invalid(G_ALLOC_WCS, xproj, message);
}
} // endif: HDU was valid
// Return
return;
}
/***********************************************************************//**
* @brief Create FITS HDU containing Healpix data
*
* This method allocates a binary table HDU that contains the Healpix data.
* Deallocation of the table has to be done by the client.
***************************************************************************/
GFitsBinTable* GSkymap::create_healpix_hdu(void) const
{
// Initialise result to NULL pointer
GFitsBinTable* hdu = NULL;
// Compute size of Healpix data
int size = m_num_pixels * m_num_maps;
// Continue only if we have pixels
if (size > 0) {
// Set number of rows and columns
int rows = m_num_pixels;
int number = m_num_maps;
// Create column to hold Healpix data
GFitsTableDoubleCol column = GFitsTableDoubleCol("DATA", rows, number);
// Fill data into column
double* ptr = m_pixels;
for (int inx = 0; inx < number; ++inx) {
for (int row = 0; row < rows; ++row) {
column(row,inx) = *ptr++;
}
}
// Create HDU that contains Healpix map in a binary table
hdu = new GFitsBinTable(rows);
hdu->append_column(column);
} // endif: there were pixels
// ... otherwise create an empty header
else {
hdu = new GFitsBinTable;
}
// Set extension name
hdu->extname("HEALPIX");
// If we have WCS information then write into FITS header
if (m_wcs != NULL) m_wcs->write(hdu);
// Set additional keywords
hdu->card("NBRBINS", m_num_maps, "Number of HEALPix maps");
// Return HDU
return hdu;
}
/***********************************************************************//**
* @brief Create FITS HDU containing WCS image
*
* This method allocates an image HDU that contains the WCS image data.
* Deallocation of the image has to be done by the client.
*
* @todo Set additional keywords.
***************************************************************************/
GFitsImageDouble* GSkymap::create_wcs_hdu(void) const
{
// Initialise result to NULL pointer
GFitsImageDouble* hdu = NULL;
// Compute size of Healpix data
int size = m_num_pixels * m_num_maps;
// Continue only if we have pixels
if (size > 0) {
// Set axis parameters for image construction
int naxis = (m_num_maps == 1) ? 2 : 3;
int naxes[] = {m_num_x, m_num_y, m_num_maps};
// Allocate image
hdu = new GFitsImageDouble(naxis, naxes);
// Store data in image
if (naxis == 2) {
double* ptr = m_pixels;
for (int iy = 0; iy < m_num_y; ++iy) {
for (int ix = 0; ix < m_num_x; ++ix) {
(*hdu)(ix,iy) = *ptr++;
}
}
}
else {
double* ptr = m_pixels;
for (int imap = 0; imap < m_num_maps; ++imap) {
for (int iy = 0; iy < m_num_y; ++iy) {
for (int ix = 0; ix < m_num_x; ++ix) {
(*hdu)(ix,iy,imap) = *ptr++;
}
}
}
}
} // endif: there were pixels
// ... otherwise create an empty header
else {
hdu = new GFitsImageDouble;
}
// Set extension name
hdu->extname("IMAGE");
// If we have WCS information then write into FITS header
if (m_wcs != NULL) m_wcs->write(hdu);
// Set additional keywords
//TODO
// Return HDU
return hdu;
}