From ac6b92b73346342caff539597a5d8b8e49db455d Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Mon, 30 May 2016 15:46:10 +0200 Subject: [PATCH] Add shape() and ndim() methods to GSkyMap. --- ChangeLog | 3 ++- NEWS | 11 ++++++++++- include/GSkyMap.hpp | 95 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------- pyext/GSkyMap.i | 71 +++++++++++++++++++++++++++++++++++++++++------------------------------ src/sky/GSkyMap.cpp | 276 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------------------------------------- test/test_GSky.cpp | 50 ++++++++++++++++++++++++++++++++++++++++++++++++++ test/test_GSky.py | 169 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------- 7 files changed, 493 insertions(+), 182 deletions(-) diff --git a/ChangeLog b/ChangeLog index aede653..44e395a 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,4 +1,4 @@ -2016-04-25 Juergen Knoedlseder +2016-05-30 Juergen Knoedlseder * Version 1.1.0 released ========================= @@ -10,6 +10,7 @@ Add GApplication::logFileClose() method. Add GCTAObservation filename constructor (#1744). Add gammalib::number() function (#1743). + Add GSkyMap::shape() and ndim() methods (#1733). Add energy dispersion cube for stacked analysis (#1721). Add GEnergy division operator. Support Fermi LAT Pass 8 response functions (#1517) diff --git a/NEWS b/NEWS index 1e8dff6..6cbff81 100644 --- a/NEWS +++ b/NEWS @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 1.1.0 Author(s) : Juergen Knoedlseder -25 April 2016 +30 May 2016 1. Introduction @@ -25,6 +25,8 @@ The following methods have been added to the public interface: - GFilename::path() - GFilename::file() - GCsv::append() +- GSkyMap::shape() +- GSkyMap::ndim() - GModelSky::spatial() - GModelSky::spectral() - GModelSky::temporal() @@ -210,6 +212,13 @@ A GSkyMap::publish() method has been added that publishes the sky map on a Virtual Observatory Hub. This finally makes GammaLib interoperable with VO imaging tools (#1231). +A shape member was added to the GSkyMap class to enable handling of +images with more than 3 dimensions. The shape() methods set or retrieve +the number of maps in each of the dimensions, the ndim() method returns +the number of dimensions. Note that a longitude and latitude dimensions +are not counted in the map dimensions, i.e. a 3D image cube with +longitude, latitude and a number of maps has a dimension of 1 (#1733). + 13. Support module ------------------ diff --git a/include/GSkyMap.hpp b/include/GSkyMap.hpp index 5328027..71f88f9 100644 --- a/include/GSkyMap.hpp +++ b/include/GSkyMap.hpp @@ -29,6 +29,7 @@ /* __ Includes ___________________________________________________________ */ #include +#include #include "GBase.hpp" #include "GSkyDir.hpp" #include "GSkyPixel.hpp" @@ -127,37 +128,44 @@ public: double operator()(const GSkyDir& dir, const int& map = 0) const; // Methods - void clear(void); - GSkyMap* clone(void) const; - std::string classname(void) const; - const int& npix(void) const; - const int& nx(void) const; - const int& ny(void) const; - const int& nmaps(void) const; - void nmaps(const int& nmaps); - GSkyPixel inx2pix(const int& index) const; - GSkyDir inx2dir(const int& index) const; - GSkyDir pix2dir(const GSkyPixel& pixel) const; - int pix2inx(const GSkyPixel& pixel) const; - int dir2inx(const GSkyDir& dir) const; - GSkyPixel dir2pix(const GSkyDir& dir) const; - double flux(const int& index, const int& map = 0) const; - double flux(const GSkyPixel& pixel, const int& map = 0) const; - double solidangle(const int& index) const; - double solidangle(const GSkyPixel& pixel) const; - bool contains(const GSkyDir& dir) const; - bool contains(const GSkyPixel& pixel) const; - const GSkyProjection* projection(void) const; - void projection(const GSkyProjection& proj); - const double* pixels(void) const; - GSkyMap extract(const int& map, const int& nmaps = 1) const; - void stack_maps(void); - void load(const GFilename& filename); - void save(const GFilename& filename, bool clobber = false) const; - void read(const GFitsHDU& hdu); - void write(GFits& file) const; - void publish(const std::string& name = "") const; - std::string print(const GChatter& chatter = NORMAL) const; + void clear(void); + GSkyMap* clone(void) const; + std::string classname(void) const; + const int& npix(void) const; + const int& nx(void) const; + const int& ny(void) const; + const int& nmaps(void) const; + void nmaps(const int& nmaps); + const std::vector& shape(void) const; + void shape(const int& s1); + void shape(const int& s1, const int& s2); + void shape(const int& s1, const int& s2, const int& s3); + void shape(const std::vector& shape); + int ndim(void) const; + GSkyPixel inx2pix(const int& index) const; + GSkyDir inx2dir(const int& index) const; + GSkyDir pix2dir(const GSkyPixel& pixel) const; + int pix2inx(const GSkyPixel& pixel) const; + int dir2inx(const GSkyDir& dir) const; + GSkyPixel dir2pix(const GSkyDir& dir) const; + double flux(const int& index, const int& map = 0) const; + double flux(const GSkyPixel& pixel, const int& map = 0) const; + double solidangle(const int& index) const; + double solidangle(const GSkyPixel& pixel) const; + bool contains(const GSkyDir& dir) const; + bool contains(const GSkyPixel& pixel) const; + const GSkyProjection* projection(void) const; + void projection(const GSkyProjection& proj); + const double* pixels(void) const; + GSkyMap extract(const int& map, const int& nmaps = 1) const; + void stack_maps(void); + void load(const GFilename& filename); + void save(const GFilename& filename, + const bool& clobber = false) const; + void read(const GFitsHDU& hdu); + void write(GFits& file) const; + void publish(const std::string& name = "") const; + std::string print(const GChatter& chatter = NORMAL) const; private: // Private methods @@ -187,6 +195,7 @@ private: int m_num_maps; //!< Number of maps (used for pixel allocation) int m_num_x; //!< Number of pixels in x direction (only 2D) int m_num_y; //!< Number of pixels in y direction (only 2D) + std::vector m_shape; //!< Shape of the maps GSkyProjection* m_proj; //!< Pointer to sky projection double* m_pixels; //!< Pointer to skymap pixels @@ -337,6 +346,30 @@ const int& GSkyMap::nmaps(void) const /***********************************************************************//** + * @brief Returns shape of maps + * + * @return Shape of maps in the sky map object. + ***************************************************************************/ +inline +const std::vector& GSkyMap::shape(void) const +{ + return m_shape; +} + + +/***********************************************************************//** + * @brief Returns dimension of maps + * + * @return Number of map dimensions. + ***************************************************************************/ +inline +int GSkyMap::ndim(void) const +{ + return m_shape.size(); +} + + +/***********************************************************************//** * @brief Returns pointer to sky projection * * @return Pointer to sky projection (NULL if no projection is defined). diff --git a/pyext/GSkyMap.i b/pyext/GSkyMap.i index ddc8f78..030f9b0 100644 --- a/pyext/GSkyMap.i +++ b/pyext/GSkyMap.i @@ -28,6 +28,10 @@ #include "GSkyMap.hpp" #include "GTools.hpp" %} +%include "std_vector.i" +namespace std { + %template(IntVector) vector; +} /***********************************************************************//** @@ -143,36 +147,43 @@ public: double operator()(const GSkyDir& dir, const int& map = 0) const; // Methods - void clear(void); - GSkyMap* clone(void) const; - std::string classname(void) const; - const int& npix(void) const; - const int& nx(void) const; - const int& ny(void) const; - const int& nmaps(void) const; - void nmaps(const int& nmaps); - GSkyPixel inx2pix(const int& index) const; - GSkyDir inx2dir(const int& index) const; - GSkyDir pix2dir(const GSkyPixel& pixel) const; - int pix2inx(const GSkyPixel& pixel) const; - int dir2inx(const GSkyDir& dir) const; - GSkyPixel dir2pix(const GSkyDir& dir) const; - double flux(const int& index, const int& map = 0) const; - double flux(const GSkyPixel& pixel, const int& map = 0) const; - double solidangle(const int& index) const; - double solidangle(const GSkyPixel& pixel) const; - bool contains(const GSkyDir& dir) const; - bool contains(const GSkyPixel& pixel) const; - const GSkyProjection* projection(void) const; - void projection(const GSkyProjection& proj); - const double* pixels(void) const; - GSkyMap extract(const int& map, const int& nmaps = 1) const; - void stack_maps(void); - void load(const GFilename& filename); - void save(const GFilename& filename, bool clobber = false) const; - void read(const GFitsHDU& hdu); - void write(GFits& file) const; - void publish(const std::string& name = "") const; + void clear(void); + GSkyMap* clone(void) const; + std::string classname(void) const; + const int& npix(void) const; + const int& nx(void) const; + const int& ny(void) const; + const int& nmaps(void) const; + void nmaps(const int& nmaps); + const std::vector& shape(void) const; + void shape(const int& s1); + void shape(const int& s1, const int& s2); + void shape(const int& s1, const int& s2, const int& s3); + void shape(const std::vector& shape); + int ndim(void) const; + GSkyPixel inx2pix(const int& index) const; + GSkyDir inx2dir(const int& index) const; + GSkyDir pix2dir(const GSkyPixel& pixel) const; + int pix2inx(const GSkyPixel& pixel) const; + int dir2inx(const GSkyDir& dir) const; + GSkyPixel dir2pix(const GSkyDir& dir) const; + double flux(const int& index, const int& map = 0) const; + double flux(const GSkyPixel& pixel, const int& map = 0) const; + double solidangle(const int& index) const; + double solidangle(const GSkyPixel& pixel) const; + bool contains(const GSkyDir& dir) const; + bool contains(const GSkyPixel& pixel) const; + const GSkyProjection* projection(void) const; + void projection(const GSkyProjection& proj); + const double* pixels(void) const; + GSkyMap extract(const int& map, const int& nmaps = 1) const; + void stack_maps(void); + void load(const GFilename& filename); + void save(const GFilename& filename, + const bool& clobber = false) const; + void read(const GFitsHDU& hdu); + void write(GFits& file) const; + void publish(const std::string& name = "") const; }; diff --git a/src/sky/GSkyMap.cpp b/src/sky/GSkyMap.cpp index 779b88f..0a2a340 100644 --- a/src/sky/GSkyMap.cpp +++ b/src/sky/GSkyMap.cpp @@ -51,6 +51,7 @@ #define G_CONSTRUCT_MAP "GSkyMap::GSkyMap(std::string&, std::string&,"\ " double&, double&, double& double&, int&, int&, int&)" #define G_NMAPS "GSkyMap::nmaps(int&)" +#define G_SHAPE "GSkyMap::shape(std::vector&)" #define G_OP_UNARY_ADD "GSkyMap::operator+=(GSkyMap&)" #define G_OP_UNARY_SUB "GSkyMap::operator-=(GSkyMap&)" #define G_OP_UNARY_MUL "GSkyMap::operator-=(GSkyMap&)" @@ -95,6 +96,8 @@ /***********************************************************************//** * @brief Void constructor + * + * Constructs an empty sky map. ***************************************************************************/ GSkyMap::GSkyMap(void) { @@ -110,6 +113,9 @@ GSkyMap::GSkyMap(void) * @brief FITS file constructor * * @param[in] filename FITS file name. + * + * Constructs a sky map by loading data from a FITS file. See the load() + * method for more information about the FITS formats that are supported. ***************************************************************************/ GSkyMap::GSkyMap(const GFilename& filename) { @@ -130,12 +136,13 @@ GSkyMap::GSkyMap(const GFilename& filename) * @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). + * @param[in] nmaps Number of maps in set. * * @exception GException::skymap_bad_par * Invalid sky map parameter. * - * Constructs sky map in Healpix pixelisation. + * Constructs @p nmaps identical all sky maps in Healpix pixelisation. All + * pixels of the sky maps will be initialised to values of zero. ***************************************************************************/ GSkyMap::GSkyMap(const std::string& coords, const int& nside, @@ -155,9 +162,10 @@ GSkyMap::GSkyMap(const std::string& coords, GHealpix* projection = new GHealpix(nside, order, coords); m_proj = projection; - // Set number of pixels and number of maps + // Set number of pixels, number of maps and shape of maps m_num_pixels = projection->npix(); m_num_maps = nmaps; + m_shape.push_back(nmaps); // Allocate pixels alloc_pixels(); @@ -183,7 +191,9 @@ GSkyMap::GSkyMap(const std::string& coords, * @exception GException::skymap_bad_par * Invalid sky map parameter. * - * Constructs sky map in World Coordinate System projection. + * Constructs @p nmaps identical all sky maps in World Coordinate System + * projection. All pixels of the sky maps will be initialised to values of + * zero. ***************************************************************************/ GSkyMap::GSkyMap(const std::string& wcs, const std::string& coords, @@ -224,11 +234,12 @@ GSkyMap::GSkyMap(const std::string& wcs, set_wcs(wcs, coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2, cd, pv2); - // Set number of pixels and number of maps + // Set number of pixels, number of maps and shape of maps m_num_x = nx; m_num_y = ny; m_num_pixels = m_num_x * m_num_y; m_num_maps = nmaps; + m_shape.push_back(nmaps); // Allocate pixels alloc_pixels(); @@ -242,6 +253,8 @@ GSkyMap::GSkyMap(const std::string& wcs, * @brief Copy constructor * * @param[in] map Sky map. + * + * Constructs sky maps by copying data from another sky map object. ***************************************************************************/ GSkyMap::GSkyMap(const GSkyMap& map) { @@ -280,6 +293,8 @@ GSkyMap::~GSkyMap(void) * * @param[in] map Sky map. * @return Sky map. + * + * Assigns one sky map to another. ***************************************************************************/ GSkyMap& GSkyMap::operator=(const GSkyMap& map) { @@ -308,7 +323,7 @@ GSkyMap& GSkyMap::operator=(const GSkyMap& map) * @param[in] value Value. * @return Sky map. * - * Sets all pixels to the specified @p value. + * Sets all sky map pixels to the specified @p value. ***************************************************************************/ GSkyMap& GSkyMap::operator=(const double& value) { @@ -335,7 +350,7 @@ GSkyMap& GSkyMap::operator=(const double& value) * Mismatch between number of maps in skymap object. * * Adds the content of @p map to the skymap. The operator only works on sky - * maps with an identical number of layers. The content is added by + * maps with an identical number of maps. The content is added by * bi-linearily interpolating the values in the source sky map, allowing thus * for a reprojection of sky map values. * @@ -1043,6 +1058,11 @@ void GSkyMap::nmaps(const int& nmaps) // Set number of maps m_num_maps = nmaps; + // Set shape of maps to linear array since we do not know how to + // arrange the new maps + m_shape.clear(); + m_shape.push_back(nmaps); + } // endif: map had pixels // Return @@ -1051,6 +1071,116 @@ void GSkyMap::nmaps(const int& nmaps) /***********************************************************************//** + * @brief Set one-dimensional shape of maps + * + * @param[in] s1 Axis length of first dimension. + ***************************************************************************/ +void GSkyMap::shape(const int& s1) +{ + // Initialise vector + std::vector shape; + + // Set vector elements + shape.push_back(s1); + + // Call vector method + this->shape(shape); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set two-dimensional shape of maps + * + * @param[in] s1 Axis length of first dimension. + * @param[in] s2 Axis length of second dimension. + ***************************************************************************/ +void GSkyMap::shape(const int& s1, const int& s2) +{ + // Initialise vector + std::vector shape; + + // Set vector elements + shape.push_back(s1); + shape.push_back(s2); + + // Call vector method + this->shape(shape); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set three-dimensional shape of maps + * + * @param[in] s1 Axis length of first dimension. + * @param[in] s2 Axis length of second dimension. + * @param[in] s3 Axis length of second dimension. + ***************************************************************************/ +void GSkyMap::shape(const int& s1, const int& s2, const int& s3) +{ + // Initialise vector + std::vector shape; + + // Set vector elements + shape.push_back(s1); + shape.push_back(s2); + shape.push_back(s3); + + // Call vector method + this->shape(shape); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set shape of maps + * + * @param[in] shape Shape vector. + * + * @exception GException::invalid_argument + * Invalid shape factorisation specified. + * + * Defines a shape for the maps in the object. The shape specifies how the + * maps are arranged in a n-dimensional array. + ***************************************************************************/ +void GSkyMap::shape(const std::vector& shape) +{ + // Computes the resulting number of maps + int nmaps = 0; + if (shape.size() > 0) { + nmaps = 1; + for (int i = 0; i < shape.size(); ++i) { + nmaps *= shape[i]; + } + } + + // Throw an exception if resulting number of maps is not equal to the + // existing number of maps + if (nmaps != m_num_maps) { + std::string msg = "Sky map shaping requires "+gammalib::str(nmaps)+ + " maps while "+gammalib::str(m_num_maps)+" maps " + "are available. Please specify an appropriate " + "factorisation or modify the number of available " + "maps."; + throw GException::invalid_argument(G_SHAPE, msg); + } + + // Set shape + m_shape = shape; + + // Return + return; +} + + +/***********************************************************************//** * @brief Converts pixel index into sky map pixel * * @param[in] index Pixel index [0,...,npix()-1]. @@ -1624,7 +1754,7 @@ bool GSkyMap::contains(const GSkyPixel& pixel) const * @brief Extract maps into a new sky map object * * @param[in] map First map to extract - * @param[in] nmaps Number of maps to extract (default: 1) + * @param[in] nmaps Number of maps to extract * @return Extracted map(s). * * @exception GException::out_of_range @@ -1685,6 +1815,10 @@ GSkyMap GSkyMap::extract(const int& map, const int& nmaps) const // Set number of maps result.m_num_maps = nmaps; + // Reset shape of maps since we do not know how to arrange the maps + result.m_shape.clear(); + result.m_shape.push_back(nmaps); + // Return map return result; } @@ -1726,6 +1860,10 @@ void GSkyMap::stack_maps(void) // Set number of maps to 1 m_num_maps = 1; + // Reset shape of maps + m_shape.clear(); + m_shape.push_back(m_num_maps); + } // endif: map had pixels // Return @@ -1839,13 +1977,13 @@ void GSkyMap::load(const GFilename& filename) * @brief Save sky map into FITS file * * @param[in] filename FITS file name. - * @param[in] clobber Overwrite existing file? (default: false) + * @param[in] clobber Overwrite existing file? * * Saves the sky map into a FITS file. If the file exists already the method * will append the sky map to the FITS file in case that the @p clobber * parameter is set to true. ***************************************************************************/ -void GSkyMap::save(const GFilename& filename, bool clobber) const +void GSkyMap::save(const GFilename& filename, const bool& clobber) const { // Continue only if we have data to save if (m_proj != NULL) { @@ -2023,11 +2161,11 @@ std::string GSkyMap::print(const GChatter& chatter) const // Append header result.append("=== GSkyMap ==="); - // Append information + // Append number of pixels and number of maps result.append("\n"+gammalib::parformat("Number of pixels")); result.append(gammalib::str(m_num_pixels)); - result.append("\n"+gammalib::parformat("Number of maps")); - result.append(gammalib::str(m_num_maps)); + + // Append WCS dimension information if (m_proj != NULL && m_proj->size() == 2) { result.append("\n"+gammalib::parformat("X axis dimension")); result.append(gammalib::str(m_num_x)); @@ -2035,6 +2173,27 @@ std::string GSkyMap::print(const GChatter& chatter) const result.append(gammalib::str(m_num_y)); } + // Append number of maps + result.append("\n"+gammalib::parformat("Number of maps")); + result.append(gammalib::str(m_num_maps)); + + // Append shape of maps + result.append("\n"+gammalib::parformat("Shape of maps")); + if (m_shape.size() > 0) { + std::string shape = "("; + for (int i = 0; i < m_shape.size(); ++i) { + if (i > 0) { + shape += ", "; + } + shape += gammalib::str(m_shape[i]); + } + shape += ")"; + result.append(shape); + } + else { + result.append("not defined"); + } + // Append sky projection information if (gammalib::reduce(chatter) > SILENT) { if (m_proj != NULL) { @@ -2068,6 +2227,7 @@ void GSkyMap::init_members(void) m_num_maps = 0; m_num_x = 0; m_num_y = 0; + m_shape.clear(); m_proj = NULL; m_pixels = NULL; @@ -2118,6 +2278,7 @@ void GSkyMap::copy_members(const GSkyMap& map) m_num_maps = map.m_num_maps; m_num_x = map.m_num_x; m_num_y = map.m_num_y; + m_shape = map.m_shape; // Copy computation cache m_hascache = map.m_hascache; @@ -2157,12 +2318,6 @@ void GSkyMap::free_members(void) m_proj = 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; } @@ -2293,6 +2448,10 @@ void GSkyMap::read_healpix(const GFitsTable& table) std::cout << "m_num_maps=" << m_num_maps << std::endl; #endif + // Initialise shape of maps + m_shape.clear(); + m_shape.push_back(m_num_maps); + // Allocate pixels to hold the map alloc_pixels(); @@ -2360,34 +2519,50 @@ void GSkyMap::read_healpix(const GFitsTable& table) * * @param[in] image FITS image. * + * @exception GException::invalid_argument + * FITS image has less than two dimensions * @exception GException::invalid_value - * WCS image has an invalid dimension or covers a too large - * range. + * Sky map covers more than 360 deg in longitude or 180 deg in + * latitude + * + * Reads sky maps from a FITS image extension containing a set of maps given + * in World Coordinate System. The method handles general n-dimensional + * images and sets the map shape attribute according to the number of map + * dimensions found in the FITS HDU. ***************************************************************************/ void GSkyMap::read_wcs(const GFitsImage& image) { + // Throw an exception if the FITS image is not at least a 2D image + if (image.naxis() < 2) { + std::string msg = "Sky map has "+gammalib::str(image.naxis())+ + " dimensions, which is less than the two dimensions" + " that are required for a WCS image."; + throw GException::invalid_argument(G_READ_WCS, msg); + } + // Allocate WCS alloc_wcs(image); // Read projection information from FITS header m_proj->read(image); - // Extract map dimension and number of maps from image + // Clear map shape + m_shape.clear(); + + // Extract map dimension, number of maps and map shape from FITS image + m_num_x = image.naxes(0); + m_num_y = image.naxes(1); if (image.naxis() == 2) { - m_num_x = image.naxes(0); - m_num_y = image.naxes(1); m_num_maps = 1; - } - else if (image.naxis() >= 3) { - m_num_x = image.naxes(0); - m_num_y = image.naxes(1); - m_num_maps = image.naxes(2); + m_shape.push_back(m_num_maps); } else { - std::string msg = "Skymap has a dimension of "+ - gammalib::str(image.naxis())+" but only dimensions " - "of 2 or 3 are supported."; - throw GException::invalid_value(G_READ_WCS, msg); + m_num_maps = image.naxes(2); + m_shape.push_back(m_num_maps); + for (int i = 3; i < image.naxis(); ++i) { + m_num_maps *= image.naxes(i); + m_shape.push_back(image.naxes(i)); + } } #if defined(G_READ_WCS_DEBUG) std::cout << "m_num_x=" << m_num_x << std::endl; @@ -2405,23 +2580,10 @@ void GSkyMap::read_wcs(const GFitsImage& image) alloc_pixels(); // Read image - if (image.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++ = image.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++ = image.pixel(ix,iy,imap); - } - } - } + double* ptr = m_pixels; + int size = m_num_pixels * m_num_maps; + for (int i = 0; i < size; ++i) { + *ptr++ = (double)image.pixel(i); } // Check if the map is too large @@ -2570,12 +2732,22 @@ GFitsImageDouble* GSkyMap::create_wcs_hdu(void) const 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}; + int naxis = (m_num_maps == 1) ? 2 : 2 + ndim(); + + // Set dimensions of all axes + int* naxes = new int[naxis]; + naxes[0] = m_num_x; + naxes[1] = m_num_y; + for (int i = 0; i < ndim(); ++i) { + naxes[2+i] = m_shape[i]; + } // Allocate image hdu = new GFitsImageDouble(naxis, naxes); + // Free dimensions of all axes + delete [] naxes; + // Store data in image if (naxis == 2) { double* ptr = m_pixels; diff --git a/test/test_GSky.cpp b/test/test_GSky.cpp index 82dd17c..0b838aa 100644 --- a/test/test_GSky.cpp +++ b/test/test_GSky.cpp @@ -842,6 +842,12 @@ void TestGSky::test_GSkyMap(void) GSkyMap map_dst("CAR", "GAL", 0.0, 0.0, -0.1, 0.1, 100, 100, 2); GSkyMap map_new("CAR", "GAL", 0.0, 0.0, -0.1, 0.1, 100, 100, 2); + // Test map dimensions and shape + test_value(map_src.nmaps(), 2, "Check that sky map contains 2 maps"); + test_value(map_src.ndim(), 1, "Check that sky map has one dimension"); + test_value(map_src.shape().size(), 1, "Check that sky map has a shape size of 1"); + test_value(map_src.shape()[0], 2, "Check that sky map has 2 maps"); + // Fill map pixels double total_src = 0.0; for (int pix = 0; pix < map_src.npix(); ++pix) { @@ -1027,6 +1033,50 @@ void TestGSky::test_GSkyMap(void) test_value(total_extract, total_src, 1.0e-3, "Test extract() method with 2 maps"); test_value(map_extract.nmaps(), 2, "Test extract() method with 2 maps"); + // Define one more map for shaping manipulation + GSkyMap map_shape("CAR", "GAL", 0.0, 0.0, -1.0, 1.0, 10, 10, 12); + + // Test initial map dimensions and shape + test_value(map_shape.nmaps(), 12, "Check that sky map contains 12 maps"); + test_value(map_shape.ndim(), 1, "Check that sky map has one dimension"); + test_value(map_shape.shape().size(), 1, "Check that sky map has a shape size of 1"); + test_value(map_shape.shape()[0], 12, "Check that sky map has 12 maps"); + + // Set new map shape + map_shape.shape(3,4); + + // Test map dimensions and shape + test_value(map_shape.nmaps(), 12, "Check that sky map contains 12 maps"); + test_value(map_shape.ndim(), 2, "Check that sky map has two dimensions"); + test_value(map_shape.shape().size(), 2, "Check that sky map has a shape size of 2"); + test_value(map_shape.shape()[0], 3, "Check that sky map has 3 maps in first dimension"); + test_value(map_shape.shape()[1], 4, "Check that sky map has 4 maps in second dimension"); + + // Set new map shape + map_shape.shape(2,3,2); + + // Test map dimensions and shape + test_value(map_shape.nmaps(), 12, "Check that sky map contains 12 maps"); + test_value(map_shape.ndim(), 3, "Check that sky map has three dimensions"); + test_value(map_shape.shape().size(), 3, "Check that sky map has a shape size of 3"); + test_value(map_shape.shape()[0], 2, "Check that sky map has 2 maps in first dimension"); + test_value(map_shape.shape()[1], 3, "Check that sky map has 3 maps in second dimension"); + test_value(map_shape.shape()[2], 2, "Check that sky map has 2 maps in third dimension"); + + // Save map + map_shape.save("test_map_shape.fits", true); + + // Load map + GSkyMap map_load_shape("test_map_shape.fits"); + + // Test map dimensions and shape + test_value(map_load_shape.nmaps(), 12, "Check that sky map contains 12 maps"); + test_value(map_load_shape.ndim(), 3, "Check that sky map has three dimensions"); + test_value(map_load_shape.shape().size(), 3, "Check that sky map has a shape size of 3"); + test_value(map_load_shape.shape()[0], 2, "Check that sky map has 2 maps in first dimension"); + test_value(map_load_shape.shape()[1], 3, "Check that sky map has 3 maps in second dimension"); + test_value(map_load_shape.shape()[2], 2, "Check that sky map has 2 maps in third dimension"); + // Exit test return; } diff --git a/test/test_GSky.py b/test/test_GSky.py index b6b5cf6..6603871 100755 --- a/test/test_GSky.py +++ b/test/test_GSky.py @@ -1,7 +1,7 @@ # ========================================================================== # This module performs unit tests for the GammaLib sky module. # -# Copyright (C) 2012-2015 Juergen Knoedlseder +# Copyright (C) 2012-2016 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 @@ -71,6 +71,8 @@ class Test(gammalib.GPythonTestSuite): self.name("sky") # Append tests + self.append(self.test_methods, "Test sky map methods") + self.append(self.test_operators, "Test sky map operators") self.append(self.test_skymap_healpix, "Test HEALPix map") self.append(self.test_skymap_ait, "Test AIT projection map") self.append(self.test_skymap_azp, "Test AZP projection map") @@ -80,11 +82,109 @@ class Test(gammalib.GPythonTestSuite): self.append(self.test_skymap_stg, "Test STG projection map") self.append(self.test_skymap_tan, "Test TAN projection map") self.append(self.test_fk5_to_galactic, "Test FK5 to Galactic coordinate conversion") - self.append(self.test_operators, "Test skymap operators") # Return return + # Test sky map methods + def test_methods(self): + """ + Test the sky map methods. + """ + # Setup sky maps + map = gammalib.GSkyMap("CAR", "CEL", 83.63, 22.01, -3.7, 2.6, 10, 8, 12) + + # Set sky map shape + map.shape(3,4) + + # Test map dimensions and shape + self.test_value(map.npix(), 80, "Check that sky map has 80 pixels"); + self.test_value(map.nx(), 10, "Check that sky map has X dimension of 10"); + self.test_value(map.ny(), 8, "Check that sky map has Y dimension of 8"); + self.test_value(map.nmaps(), 12, "Check that sky map contains 12 maps"); + self.test_value(map.ndim(), 2, "Check that sky map has two dimensions"); + self.test_value(len(map.shape()), 2, "Check that sky map has a shape size of 2"); + self.test_value(map.shape()[0], 3, "Check that sky map has 3 maps in first dimension"); + self.test_value(map.shape()[1], 4, "Check that sky map has 4 maps in second dimension"); + + # Set new map shape + map.shape(2,3,2); + + # Test map dimensions and shape + self.test_value(map.nmaps(), 12, "Check that sky map contains 12 maps"); + self.test_value(map.ndim(), 3, "Check that sky map has three dimensions"); + self.test_value(len(map.shape()), 3, "Check that sky map has a shape size of 3"); + self.test_value(map.shape()[0], 2, "Check that sky map has 2 maps in first dimension"); + self.test_value(map.shape()[1], 3, "Check that sky map has 3 maps in second dimension"); + self.test_value(map.shape()[2], 2, "Check that sky map has 2 maps in third dimension"); + + # Return + return + + # Test sky map operators + def test_operators(self): + """ + Test the sky map operators. + """ + # Setup sky maps + map = gammalib.GSkyMap("CAR", "CEL", 83.6331, 22.0145, -3.7, 2.6, 2, 2) + map[0] = 1.0 + map[1] = 2.0 + map[2] = 3.0 + map[3] = 4.0 + map_b = map.copy() + + # Addition operator + map += map_b + self.test_value(map[0], 2.0) + self.test_value(map[1], 4.0) + self.test_value(map[2], 6.0) + self.test_value(map[3], 8.0) + + # Multiplication operator + map *= map_b + self.test_value(map[0], 2.0) + self.test_value(map[1], 8.0) + self.test_value(map[2], 18.0) + self.test_value(map[3], 32.0) + + # Subtraction operator + map -= map_b + self.test_value(map[0], 1.0) + self.test_value(map[1], 6.0) + self.test_value(map[2], 15.0) + self.test_value(map[3], 28.0) + + # Division operator + map /= map_b + self.test_value(map[0], 1.0) + self.test_value(map[1], 3.0) + self.test_value(map[2], 5.0) + self.test_value(map[3], 7.0) + + # Scaling operator + map *= 2.0 + self.test_value(map[0], 2.0) + self.test_value(map[1], 6.0) + self.test_value(map[2], 10.0) + self.test_value(map[3], 14.0) + + # Division operator + map /= 2.0 + self.test_value(map[0], 1.0) + self.test_value(map[1], 3.0) + self.test_value(map[2], 5.0) + self.test_value(map[3], 7.0) + + # Access operator (tests also proper iteration) + sum = 0.0 + for pix in map: + sum += pix + self.test_value(sum, 16.0) + + # Return + return + # Generic skymap pixel transformation test def test_skymap_pixels(self, pixels, map): """ @@ -315,68 +415,3 @@ class Test(gammalib.GPythonTestSuite): # Return return - - # Test skymap operators - def test_operators(self): - """ - Test the skymap operators. - """ - # Setup skymaps - map = gammalib.GSkyMap("CAR", "CEL", 83.6331, 22.0145, -3.7, 2.6, 2, 2) - map[0] = 1.0 - map[1] = 2.0 - map[2] = 3.0 - map[3] = 4.0 - map_b = map.copy() - - # Addition operator - map += map_b - self.test_value(map[0], 2.0) - self.test_value(map[1], 4.0) - self.test_value(map[2], 6.0) - self.test_value(map[3], 8.0) - - # Multiplication operator - map *= map_b - self.test_value(map[0], 2.0) - self.test_value(map[1], 8.0) - self.test_value(map[2], 18.0) - self.test_value(map[3], 32.0) - - # Subtraction operator - map -= map_b - self.test_value(map[0], 1.0) - self.test_value(map[1], 6.0) - self.test_value(map[2], 15.0) - self.test_value(map[3], 28.0) - - # Division operator - map /= map_b - self.test_value(map[0], 1.0) - self.test_value(map[1], 3.0) - self.test_value(map[2], 5.0) - self.test_value(map[3], 7.0) - - # Scaling operator - map *= 2.0 - self.test_value(map[0], 2.0) - self.test_value(map[1], 6.0) - self.test_value(map[2], 10.0) - self.test_value(map[3], 14.0) - - # Division operator - map /= 2.0 - self.test_value(map[0], 1.0) - self.test_value(map[1], 3.0) - self.test_value(map[2], 5.0) - self.test_value(map[3], 7.0) - - # Access operator (tests also proper iteration) - sum = 0.0 - for pix in map: - sum += pix - self.test_value(sum, 16.0) - - # Return - return - \ No newline at end of file -- libgit2 0.21.4