Eirik Hatlen / gammalib

Commit ac6b92b73346342caff539597a5d8b8e49db455d

Authored by Jürgen Knödlseder
1 parent 9f845a0c

Add shape() and ndim() methods to GSkyMap. 

The shape() and ndim() methods allow the handling of images with more than three dimensions in GSkyMap. Now higher dimensional images are properly read and written by the class, and the shape() method allows setting of higher dimensional map parameters.
Inline Side-by-side
Showing 7 changed files with 493 additions and 182 deletions
ChangeLog
  View file @ac6b92b
1   -2016-04-25 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
  1 +2016-05-30 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
2 2  
3 3 * Version 1.1.0 released
4 4 =========================
... ... @@ -10,6 +10,7 @@
10 10 Add GApplication::logFileClose() method.
11 11 Add GCTAObservation filename constructor (#1744).
12 12 Add gammalib::number() function (#1743).
  13 + Add GSkyMap::shape() and ndim() methods (#1733).
13 14 Add energy dispersion cube for stacked analysis (#1721).
14 15 Add GEnergy division operator.
15 16 Support Fermi LAT Pass 8 response functions (#1517)
... ...
... ... @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 1.1.0
2 2  
3 3 Author(s) : Juergen Knoedlseder
4 4  
5   -25 April 2016
  5 +30 May 2016
6 6  
7 7  
8 8 1. Introduction
... ... @@ -25,6 +25,8 @@ The following methods have been added to the public interface:
25 25 - GFilename::path()
26 26 - GFilename::file()
27 27 - GCsv::append()
  28 +- GSkyMap::shape()
  29 +- GSkyMap::ndim()
28 30 - GModelSky::spatial()
29 31 - GModelSky::spectral()
30 32 - GModelSky::temporal()
... ... @@ -210,6 +212,13 @@ A GSkyMap::publish() method has been added that publishes the sky map
210 212 on a Virtual Observatory Hub. This finally makes GammaLib interoperable
211 213 with VO imaging tools (#1231).
212 214  
  215 +A shape member was added to the GSkyMap class to enable handling of
  216 +images with more than 3 dimensions. The shape() methods set or retrieve
  217 +the number of maps in each of the dimensions, the ndim() method returns
  218 +the number of dimensions. Note that a longitude and latitude dimensions
  219 +are not counted in the map dimensions, i.e. a 3D image cube with
  220 +longitude, latitude and a number of maps has a dimension of 1 (#1733).
  221 +
213 222  
214 223 13. Support module
215 224 ------------------
... ...
include/GSkyMap.hpp
  View file @ac6b92b
... ... @@ -29,6 +29,7 @@
29 29  
30 30 /* __ Includes ___________________________________________________________ */
31 31 #include <string>
  32 +#include <vector>
32 33 #include "GBase.hpp"
33 34 #include "GSkyDir.hpp"
34 35 #include "GSkyPixel.hpp"
... ... @@ -127,37 +128,44 @@ public:
127 128 double operator()(const GSkyDir& dir, const int& map = 0) const;
128 129  
129 130 // Methods
130   - void clear(void);
131   - GSkyMap* clone(void) const;
132   - std::string classname(void) const;
133   - const int& npix(void) const;
134   - const int& nx(void) const;
135   - const int& ny(void) const;
136   - const int& nmaps(void) const;
137   - void nmaps(const int& nmaps);
138   - GSkyPixel inx2pix(const int& index) const;
139   - GSkyDir inx2dir(const int& index) const;
140   - GSkyDir pix2dir(const GSkyPixel& pixel) const;
141   - int pix2inx(const GSkyPixel& pixel) const;
142   - int dir2inx(const GSkyDir& dir) const;
143   - GSkyPixel dir2pix(const GSkyDir& dir) const;
144   - double flux(const int& index, const int& map = 0) const;
145   - double flux(const GSkyPixel& pixel, const int& map = 0) const;
146   - double solidangle(const int& index) const;
147   - double solidangle(const GSkyPixel& pixel) const;
148   - bool contains(const GSkyDir& dir) const;
149   - bool contains(const GSkyPixel& pixel) const;
150   - const GSkyProjection* projection(void) const;
151   - void projection(const GSkyProjection& proj);
152   - const double* pixels(void) const;
153   - GSkyMap extract(const int& map, const int& nmaps = 1) const;
154   - void stack_maps(void);
155   - void load(const GFilename& filename);
156   - void save(const GFilename& filename, bool clobber = false) const;
157   - void read(const GFitsHDU& hdu);
158   - void write(GFits& file) const;
159   - void publish(const std::string& name = "") const;
160   - std::string print(const GChatter& chatter = NORMAL) const;
  131 + void clear(void);
  132 + GSkyMap* clone(void) const;
  133 + std::string classname(void) const;
  134 + const int& npix(void) const;
  135 + const int& nx(void) const;
  136 + const int& ny(void) const;
  137 + const int& nmaps(void) const;
  138 + void nmaps(const int& nmaps);
  139 + const std::vector<int>& shape(void) const;
  140 + void shape(const int& s1);
  141 + void shape(const int& s1, const int& s2);
  142 + void shape(const int& s1, const int& s2, const int& s3);
  143 + void shape(const std::vector<int>& shape);
  144 + int ndim(void) const;
  145 + GSkyPixel inx2pix(const int& index) const;
  146 + GSkyDir inx2dir(const int& index) const;
  147 + GSkyDir pix2dir(const GSkyPixel& pixel) const;
  148 + int pix2inx(const GSkyPixel& pixel) const;
  149 + int dir2inx(const GSkyDir& dir) const;
  150 + GSkyPixel dir2pix(const GSkyDir& dir) const;
  151 + double flux(const int& index, const int& map = 0) const;
  152 + double flux(const GSkyPixel& pixel, const int& map = 0) const;
  153 + double solidangle(const int& index) const;
  154 + double solidangle(const GSkyPixel& pixel) const;
  155 + bool contains(const GSkyDir& dir) const;
  156 + bool contains(const GSkyPixel& pixel) const;
  157 + const GSkyProjection* projection(void) const;
  158 + void projection(const GSkyProjection& proj);
  159 + const double* pixels(void) const;
  160 + GSkyMap extract(const int& map, const int& nmaps = 1) const;
  161 + void stack_maps(void);
  162 + void load(const GFilename& filename);
  163 + void save(const GFilename& filename,
  164 + const bool& clobber = false) const;
  165 + void read(const GFitsHDU& hdu);
  166 + void write(GFits& file) const;
  167 + void publish(const std::string& name = "") const;
  168 + std::string print(const GChatter& chatter = NORMAL) const;
161 169  
162 170 private:
163 171 // Private methods
... ... @@ -187,6 +195,7 @@ private:
187 195 int m_num_maps; //!< Number of maps (used for pixel allocation)
188 196 int m_num_x; //!< Number of pixels in x direction (only 2D)
189 197 int m_num_y; //!< Number of pixels in y direction (only 2D)
  198 + std::vector<int> m_shape; //!< Shape of the maps
190 199 GSkyProjection* m_proj; //!< Pointer to sky projection
191 200 double* m_pixels; //!< Pointer to skymap pixels
192 201  
... ... @@ -337,6 +346,30 @@ const int&amp; GSkyMap::nmaps(void) const
337 346  
338 347  
339 348 /***********************************************************************//**
  349 + * @brief Returns shape of maps
  350 + *
  351 + * @return Shape of maps in the sky map object.
  352 + ***************************************************************************/
  353 +inline
  354 +const std::vector<int>& GSkyMap::shape(void) const
  355 +{
  356 + return m_shape;
  357 +}
  358 +
  359 +
  360 +/***********************************************************************//**
  361 + * @brief Returns dimension of maps
  362 + *
  363 + * @return Number of map dimensions.
  364 + ***************************************************************************/
  365 +inline
  366 +int GSkyMap::ndim(void) const
  367 +{
  368 + return m_shape.size();
  369 +}
  370 +
  371 +
  372 +/***********************************************************************//**
340 373 * @brief Returns pointer to sky projection
341 374 *
342 375 * @return Pointer to sky projection (NULL if no projection is defined).
... ...
pyext/GSkyMap.i
  View file @ac6b92b
... ... @@ -28,6 +28,10 @@
28 28 #include "GSkyMap.hpp"
29 29 #include "GTools.hpp"
30 30 %}
  31 +%include "std_vector.i"
  32 +namespace std {
  33 + %template(IntVector) vector<int>;
  34 +}
31 35  
32 36  
33 37 /***********************************************************************//**
... ... @@ -143,36 +147,43 @@ public:
143 147 double operator()(const GSkyDir& dir, const int& map = 0) const;
144 148  
145 149 // Methods
146   - void clear(void);
147   - GSkyMap* clone(void) const;
148   - std::string classname(void) const;
149   - const int& npix(void) const;
150   - const int& nx(void) const;
151   - const int& ny(void) const;
152   - const int& nmaps(void) const;
153   - void nmaps(const int& nmaps);
154   - GSkyPixel inx2pix(const int& index) const;
155   - GSkyDir inx2dir(const int& index) const;
156   - GSkyDir pix2dir(const GSkyPixel& pixel) const;
157   - int pix2inx(const GSkyPixel& pixel) const;
158   - int dir2inx(const GSkyDir& dir) const;
159   - GSkyPixel dir2pix(const GSkyDir& dir) const;
160   - double flux(const int& index, const int& map = 0) const;
161   - double flux(const GSkyPixel& pixel, const int& map = 0) const;
162   - double solidangle(const int& index) const;
163   - double solidangle(const GSkyPixel& pixel) const;
164   - bool contains(const GSkyDir& dir) const;
165   - bool contains(const GSkyPixel& pixel) const;
166   - const GSkyProjection* projection(void) const;
167   - void projection(const GSkyProjection& proj);
168   - const double* pixels(void) const;
169   - GSkyMap extract(const int& map, const int& nmaps = 1) const;
170   - void stack_maps(void);
171   - void load(const GFilename& filename);
172   - void save(const GFilename& filename, bool clobber = false) const;
173   - void read(const GFitsHDU& hdu);
174   - void write(GFits& file) const;
175   - void publish(const std::string& name = "") const;
  150 + void clear(void);
  151 + GSkyMap* clone(void) const;
  152 + std::string classname(void) const;
  153 + const int& npix(void) const;
  154 + const int& nx(void) const;
  155 + const int& ny(void) const;
  156 + const int& nmaps(void) const;
  157 + void nmaps(const int& nmaps);
  158 + const std::vector<int>& shape(void) const;
  159 + void shape(const int& s1);
  160 + void shape(const int& s1, const int& s2);
  161 + void shape(const int& s1, const int& s2, const int& s3);
  162 + void shape(const std::vector<int>& shape);
  163 + int ndim(void) const;
  164 + GSkyPixel inx2pix(const int& index) const;
  165 + GSkyDir inx2dir(const int& index) const;
  166 + GSkyDir pix2dir(const GSkyPixel& pixel) const;
  167 + int pix2inx(const GSkyPixel& pixel) const;
  168 + int dir2inx(const GSkyDir& dir) const;
  169 + GSkyPixel dir2pix(const GSkyDir& dir) const;
  170 + double flux(const int& index, const int& map = 0) const;
  171 + double flux(const GSkyPixel& pixel, const int& map = 0) const;
  172 + double solidangle(const int& index) const;
  173 + double solidangle(const GSkyPixel& pixel) const;
  174 + bool contains(const GSkyDir& dir) const;
  175 + bool contains(const GSkyPixel& pixel) const;
  176 + const GSkyProjection* projection(void) const;
  177 + void projection(const GSkyProjection& proj);
  178 + const double* pixels(void) const;
  179 + GSkyMap extract(const int& map, const int& nmaps = 1) const;
  180 + void stack_maps(void);
  181 + void load(const GFilename& filename);
  182 + void save(const GFilename& filename,
  183 + const bool& clobber = false) const;
  184 + void read(const GFitsHDU& hdu);
  185 + void write(GFits& file) const;
  186 + void publish(const std::string& name = "") const;
176 187 };
177 188  
178 189  
... ...
src/sky/GSkyMap.cpp
  View file @ac6b92b
... ... @@ -51,6 +51,7 @@
51 51 #define G_CONSTRUCT_MAP "GSkyMap::GSkyMap(std::string&, std::string&,"\
52 52 " double&, double&, double& double&, int&, int&, int&)"
53 53 #define G_NMAPS "GSkyMap::nmaps(int&)"
  54 +#define G_SHAPE "GSkyMap::shape(std::vector<int>&)"
54 55 #define G_OP_UNARY_ADD "GSkyMap::operator+=(GSkyMap&)"
55 56 #define G_OP_UNARY_SUB "GSkyMap::operator-=(GSkyMap&)"
56 57 #define G_OP_UNARY_MUL "GSkyMap::operator-=(GSkyMap&)"
... ... @@ -95,6 +96,8 @@
95 96  
96 97 /***********************************************************************//**
97 98 * @brief Void constructor
  99 + *
  100 + * Constructs an empty sky map.
98 101 ***************************************************************************/
99 102 GSkyMap::GSkyMap(void)
100 103 {
... ... @@ -110,6 +113,9 @@ GSkyMap::GSkyMap(void)
110 113 * @brief FITS file constructor
111 114 *
112 115 * @param[in] filename FITS file name.
  116 + *
  117 + * Constructs a sky map by loading data from a FITS file. See the load()
  118 + * method for more information about the FITS formats that are supported.
113 119 ***************************************************************************/
114 120 GSkyMap::GSkyMap(const GFilename& filename)
115 121 {
... ... @@ -130,12 +136,13 @@ GSkyMap::GSkyMap(const GFilename&amp; filename)
130 136 * @param[in] coords Coordinate System (CEL or GAL).
131 137 * @param[in] nside Nside parameter.
132 138 * @param[in] order Pixel ordering (RING or NEST).
133   - * @param[in] nmaps Number of maps in set (default=1).
  139 + * @param[in] nmaps Number of maps in set.
134 140 *
135 141 * @exception GException::skymap_bad_par
136 142 * Invalid sky map parameter.
137 143 *
138   - * Constructs sky map in Healpix pixelisation.
  144 + * Constructs @p nmaps identical all sky maps in Healpix pixelisation. All
  145 + * pixels of the sky maps will be initialised to values of zero.
139 146 ***************************************************************************/
140 147 GSkyMap::GSkyMap(const std::string& coords,
141 148 const int& nside,
... ... @@ -155,9 +162,10 @@ GSkyMap::GSkyMap(const std::string&amp; coords,
155 162 GHealpix* projection = new GHealpix(nside, order, coords);
156 163 m_proj = projection;
157 164  
158   - // Set number of pixels and number of maps
  165 + // Set number of pixels, number of maps and shape of maps
159 166 m_num_pixels = projection->npix();
160 167 m_num_maps = nmaps;
  168 + m_shape.push_back(nmaps);
161 169  
162 170 // Allocate pixels
163 171 alloc_pixels();
... ... @@ -183,7 +191,9 @@ GSkyMap::GSkyMap(const std::string&amp; coords,
183 191 * @exception GException::skymap_bad_par
184 192 * Invalid sky map parameter.
185 193 *
186   - * Constructs sky map in World Coordinate System projection.
  194 + * Constructs @p nmaps identical all sky maps in World Coordinate System
  195 + * projection. All pixels of the sky maps will be initialised to values of
  196 + * zero.
187 197 ***************************************************************************/
188 198 GSkyMap::GSkyMap(const std::string& wcs,
189 199 const std::string& coords,
... ... @@ -224,11 +234,12 @@ GSkyMap::GSkyMap(const std::string&amp; wcs,
224 234 set_wcs(wcs, coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2,
225 235 cd, pv2);
226 236  
227   - // Set number of pixels and number of maps
  237 + // Set number of pixels, number of maps and shape of maps
228 238 m_num_x = nx;
229 239 m_num_y = ny;
230 240 m_num_pixels = m_num_x * m_num_y;
231 241 m_num_maps = nmaps;
  242 + m_shape.push_back(nmaps);
232 243  
233 244 // Allocate pixels
234 245 alloc_pixels();
... ... @@ -242,6 +253,8 @@ GSkyMap::GSkyMap(const std::string&amp; wcs,
242 253 * @brief Copy constructor
243 254 *
244 255 * @param[in] map Sky map.
  256 + *
  257 + * Constructs sky maps by copying data from another sky map object.
245 258 ***************************************************************************/
246 259 GSkyMap::GSkyMap(const GSkyMap& map)
247 260 {
... ... @@ -280,6 +293,8 @@ GSkyMap::~GSkyMap(void)
280 293 *
281 294 * @param[in] map Sky map.
282 295 * @return Sky map.
  296 + *
  297 + * Assigns one sky map to another.
283 298 ***************************************************************************/
284 299 GSkyMap& GSkyMap::operator=(const GSkyMap& map)
285 300 {
... ... @@ -308,7 +323,7 @@ GSkyMap&amp; GSkyMap::operator=(const GSkyMap&amp; map)
308 323 * @param[in] value Value.
309 324 * @return Sky map.
310 325 *
311   - * Sets all pixels to the specified @p value.
  326 + * Sets all sky map pixels to the specified @p value.
312 327 ***************************************************************************/
313 328 GSkyMap& GSkyMap::operator=(const double& value)
314 329 {
... ... @@ -335,7 +350,7 @@ GSkyMap&amp; GSkyMap::operator=(const double&amp; value)
335 350 * Mismatch between number of maps in skymap object.
336 351 *
337 352 * Adds the content of @p map to the skymap. The operator only works on sky
338   - * maps with an identical number of layers. The content is added by
  353 + * maps with an identical number of maps. The content is added by
339 354 * bi-linearily interpolating the values in the source sky map, allowing thus
340 355 * for a reprojection of sky map values.
341 356 *
... ... @@ -1043,6 +1058,11 @@ void GSkyMap::nmaps(const int&amp; nmaps)
1043 1058 // Set number of maps
1044 1059 m_num_maps = nmaps;
1045 1060  
  1061 + // Set shape of maps to linear array since we do not know how to
  1062 + // arrange the new maps
  1063 + m_shape.clear();
  1064 + m_shape.push_back(nmaps);
  1065 +
1046 1066 } // endif: map had pixels
1047 1067  
1048 1068 // Return
... ... @@ -1051,6 +1071,116 @@ void GSkyMap::nmaps(const int&amp; nmaps)
1051 1071  
1052 1072  
1053 1073 /***********************************************************************//**
  1074 + * @brief Set one-dimensional shape of maps
  1075 + *
  1076 + * @param[in] s1 Axis length of first dimension.
  1077 + ***************************************************************************/
  1078 +void GSkyMap::shape(const int& s1)
  1079 +{
  1080 + // Initialise vector
  1081 + std::vector<int> shape;
  1082 +
  1083 + // Set vector elements
  1084 + shape.push_back(s1);
  1085 +
  1086 + // Call vector method
  1087 + this->shape(shape);
  1088 +
  1089 + // Return
  1090 + return;
  1091 +}
  1092 +
  1093 +
  1094 +/***********************************************************************//**
  1095 + * @brief Set two-dimensional shape of maps
  1096 + *
  1097 + * @param[in] s1 Axis length of first dimension.
  1098 + * @param[in] s2 Axis length of second dimension.
  1099 + ***************************************************************************/
  1100 +void GSkyMap::shape(const int& s1, const int& s2)
  1101 +{
  1102 + // Initialise vector
  1103 + std::vector<int> shape;
  1104 +
  1105 + // Set vector elements
  1106 + shape.push_back(s1);
  1107 + shape.push_back(s2);
  1108 +
  1109 + // Call vector method
  1110 + this->shape(shape);
  1111 +
  1112 + // Return
  1113 + return;
  1114 +}
  1115 +
  1116 +
  1117 +/***********************************************************************//**
  1118 + * @brief Set three-dimensional shape of maps
  1119 + *
  1120 + * @param[in] s1 Axis length of first dimension.
  1121 + * @param[in] s2 Axis length of second dimension.
  1122 + * @param[in] s3 Axis length of second dimension.
  1123 + ***************************************************************************/
  1124 +void GSkyMap::shape(const int& s1, const int& s2, const int& s3)
  1125 +{
  1126 + // Initialise vector
  1127 + std::vector<int> shape;
  1128 +
  1129 + // Set vector elements
  1130 + shape.push_back(s1);
  1131 + shape.push_back(s2);
  1132 + shape.push_back(s3);
  1133 +
  1134 + // Call vector method
  1135 + this->shape(shape);
  1136 +
  1137 + // Return
  1138 + return;
  1139 +}
  1140 +
  1141 +
  1142 +/***********************************************************************//**
  1143 + * @brief Set shape of maps
  1144 + *
  1145 + * @param[in] shape Shape vector.
  1146 + *
  1147 + * @exception GException::invalid_argument
  1148 + * Invalid shape factorisation specified.
  1149 + *
  1150 + * Defines a shape for the maps in the object. The shape specifies how the
  1151 + * maps are arranged in a n-dimensional array.
  1152 + ***************************************************************************/
  1153 +void GSkyMap::shape(const std::vector<int>& shape)
  1154 +{
  1155 + // Computes the resulting number of maps
  1156 + int nmaps = 0;
  1157 + if (shape.size() > 0) {
  1158 + nmaps = 1;
  1159 + for (int i = 0; i < shape.size(); ++i) {
  1160 + nmaps *= shape[i];
  1161 + }
  1162 + }
  1163 +
  1164 + // Throw an exception if resulting number of maps is not equal to the
  1165 + // existing number of maps
  1166 + if (nmaps != m_num_maps) {
  1167 + std::string msg = "Sky map shaping requires "+gammalib::str(nmaps)+
  1168 + " maps while "+gammalib::str(m_num_maps)+" maps "
  1169 + "are available. Please specify an appropriate "
  1170 + "factorisation or modify the number of available "
  1171 + "maps.";
  1172 + throw GException::invalid_argument(G_SHAPE, msg);
  1173 + }
  1174 +
  1175 + // Set shape
  1176 + m_shape = shape;
  1177 +
  1178 + // Return
  1179 + return;
  1180 +}
  1181 +
  1182 +
  1183 +/***********************************************************************//**
1054 1184 * @brief Converts pixel index into sky map pixel
1055 1185 *
1056 1186 * @param[in] index Pixel index [0,...,npix()-1].
... ... @@ -1624,7 +1754,7 @@ bool GSkyMap::contains(const GSkyPixel&amp; pixel) const
1624 1754 * @brief Extract maps into a new sky map object
1625 1755 *
1626 1756 * @param[in] map First map to extract
1627   - * @param[in] nmaps Number of maps to extract (default: 1)
  1757 + * @param[in] nmaps Number of maps to extract
1628 1758 * @return Extracted map(s).
1629 1759 *
1630 1760 * @exception GException::out_of_range
... ... @@ -1685,6 +1815,10 @@ GSkyMap GSkyMap::extract(const int&amp; map, const int&amp; nmaps) const
1685 1815 // Set number of maps
1686 1816 result.m_num_maps = nmaps;
1687 1817  
  1818 + // Reset shape of maps since we do not know how to arrange the maps
  1819 + result.m_shape.clear();
  1820 + result.m_shape.push_back(nmaps);
  1821 +
1688 1822 // Return map
1689 1823 return result;
1690 1824 }
... ... @@ -1726,6 +1860,10 @@ void GSkyMap::stack_maps(void)
1726 1860 // Set number of maps to 1
1727 1861 m_num_maps = 1;
1728 1862  
  1863 + // Reset shape of maps
  1864 + m_shape.clear();
  1865 + m_shape.push_back(m_num_maps);
  1866 +
1729 1867 } // endif: map had pixels
1730 1868  
1731 1869 // Return
... ... @@ -1839,13 +1977,13 @@ void GSkyMap::load(const GFilename&amp; filename)
1839 1977 * @brief Save sky map into FITS file
1840 1978 *
1841 1979 * @param[in] filename FITS file name.
1842   - * @param[in] clobber Overwrite existing file? (default: false)
  1980 + * @param[in] clobber Overwrite existing file?
1843 1981 *
1844 1982 * Saves the sky map into a FITS file. If the file exists already the method
1845 1983 * will append the sky map to the FITS file in case that the @p clobber
1846 1984 * parameter is set to true.
1847 1985 ***************************************************************************/
1848   -void GSkyMap::save(const GFilename& filename, bool clobber) const
  1986 +void GSkyMap::save(const GFilename& filename, const bool& clobber) const
1849 1987 {
1850 1988 // Continue only if we have data to save
1851 1989 if (m_proj != NULL) {
... ... @@ -2023,11 +2161,11 @@ std::string GSkyMap::print(const GChatter&amp; chatter) const
2023 2161 // Append header
2024 2162 result.append("=== GSkyMap ===");
2025 2163  
2026   - // Append information
  2164 + // Append number of pixels and number of maps
2027 2165 result.append("\n"+gammalib::parformat("Number of pixels"));
2028 2166 result.append(gammalib::str(m_num_pixels));
2029   - result.append("\n"+gammalib::parformat("Number of maps"));
2030   - result.append(gammalib::str(m_num_maps));
  2167 +
  2168 + // Append WCS dimension information
2031 2169 if (m_proj != NULL && m_proj->size() == 2) {
2032 2170 result.append("\n"+gammalib::parformat("X axis dimension"));
2033 2171 result.append(gammalib::str(m_num_x));
... ... @@ -2035,6 +2173,27 @@ std::string GSkyMap::print(const GChatter&amp; chatter) const
2035 2173 result.append(gammalib::str(m_num_y));
2036 2174 }
2037 2175  
  2176 + // Append number of maps
  2177 + result.append("\n"+gammalib::parformat("Number of maps"));
  2178 + result.append(gammalib::str(m_num_maps));
  2179 +
  2180 + // Append shape of maps
  2181 + result.append("\n"+gammalib::parformat("Shape of maps"));
  2182 + if (m_shape.size() > 0) {
  2183 + std::string shape = "(";
  2184 + for (int i = 0; i < m_shape.size(); ++i) {
  2185 + if (i > 0) {
  2186 + shape += ", ";
  2187 + }
  2188 + shape += gammalib::str(m_shape[i]);
  2189 + }
  2190 + shape += ")";
  2191 + result.append(shape);
  2192 + }
  2193 + else {
  2194 + result.append("not defined");
  2195 + }
  2196 +
2038 2197 // Append sky projection information
2039 2198 if (gammalib::reduce(chatter) > SILENT) {
2040 2199 if (m_proj != NULL) {
... ... @@ -2068,6 +2227,7 @@ void GSkyMap::init_members(void)
2068 2227 m_num_maps = 0;
2069 2228 m_num_x = 0;
2070 2229 m_num_y = 0;
  2230 + m_shape.clear();
2071 2231 m_proj = NULL;
2072 2232 m_pixels = NULL;
2073 2233  
... ... @@ -2118,6 +2278,7 @@ void GSkyMap::copy_members(const GSkyMap&amp; map)
2118 2278 m_num_maps = map.m_num_maps;
2119 2279 m_num_x = map.m_num_x;
2120 2280 m_num_y = map.m_num_y;
  2281 + m_shape = map.m_shape;
2121 2282  
2122 2283 // Copy computation cache
2123 2284 m_hascache = map.m_hascache;
... ... @@ -2157,12 +2318,6 @@ void GSkyMap::free_members(void)
2157 2318 m_proj = NULL;
2158 2319 m_pixels = NULL;
2159 2320  
2160   - // Reset number of pixels
2161   - m_num_pixels = 0;
2162   - m_num_maps = 0;
2163   - m_num_x = 0;
2164   - m_num_y = 0;
2165   -
2166 2321 // Return
2167 2322 return;
2168 2323 }
... ... @@ -2293,6 +2448,10 @@ void GSkyMap::read_healpix(const GFitsTable&amp; table)
2293 2448 std::cout << "m_num_maps=" << m_num_maps << std::endl;
2294 2449 #endif
2295 2450  
  2451 + // Initialise shape of maps
  2452 + m_shape.clear();
  2453 + m_shape.push_back(m_num_maps);
  2454 +
2296 2455 // Allocate pixels to hold the map
2297 2456 alloc_pixels();
2298 2457  
... ... @@ -2360,34 +2519,50 @@ void GSkyMap::read_healpix(const GFitsTable&amp; table)
2360 2519 *
2361 2520 * @param[in] image FITS image.
2362 2521 *
  2522 + * @exception GException::invalid_argument
  2523 + * FITS image has less than two dimensions
2363 2524 * @exception GException::invalid_value
2364   - * WCS image has an invalid dimension or covers a too large
2365   - * range.
  2525 + * Sky map covers more than 360 deg in longitude or 180 deg in
  2526 + * latitude
  2527 + *
  2528 + * Reads sky maps from a FITS image extension containing a set of maps given
  2529 + * in World Coordinate System. The method handles general n-dimensional
  2530 + * images and sets the map shape attribute according to the number of map
  2531 + * dimensions found in the FITS HDU.
2366 2532 ***************************************************************************/
2367 2533 void GSkyMap::read_wcs(const GFitsImage& image)
2368 2534 {
  2535 + // Throw an exception if the FITS image is not at least a 2D image
  2536 + if (image.naxis() < 2) {
  2537 + std::string msg = "Sky map has "+gammalib::str(image.naxis())+
  2538 + " dimensions, which is less than the two dimensions"
  2539 + " that are required for a WCS image.";
  2540 + throw GException::invalid_argument(G_READ_WCS, msg);
  2541 + }
  2542 +
2369 2543 // Allocate WCS
2370 2544 alloc_wcs(image);
2371 2545  
2372 2546 // Read projection information from FITS header
2373 2547 m_proj->read(image);
2374 2548  
2375   - // Extract map dimension and number of maps from image
  2549 + // Clear map shape
  2550 + m_shape.clear();
  2551 +
  2552 + // Extract map dimension, number of maps and map shape from FITS image
  2553 + m_num_x = image.naxes(0);
  2554 + m_num_y = image.naxes(1);
2376 2555 if (image.naxis() == 2) {
2377   - m_num_x = image.naxes(0);
2378   - m_num_y = image.naxes(1);
2379 2556 m_num_maps = 1;
2380   - }
2381   - else if (image.naxis() >= 3) {
2382   - m_num_x = image.naxes(0);
2383   - m_num_y = image.naxes(1);
2384   - m_num_maps = image.naxes(2);
  2557 + m_shape.push_back(m_num_maps);
2385 2558 }
2386 2559 else {
2387   - std::string msg = "Skymap has a dimension of "+
2388   - gammalib::str(image.naxis())+" but only dimensions "
2389   - "of 2 or 3 are supported.";
2390   - throw GException::invalid_value(G_READ_WCS, msg);
  2560 + m_num_maps = image.naxes(2);
  2561 + m_shape.push_back(m_num_maps);
  2562 + for (int i = 3; i < image.naxis(); ++i) {
  2563 + m_num_maps *= image.naxes(i);
  2564 + m_shape.push_back(image.naxes(i));
  2565 + }
2391 2566 }
2392 2567 #if defined(G_READ_WCS_DEBUG)
2393 2568 std::cout << "m_num_x=" << m_num_x << std::endl;
... ... @@ -2405,23 +2580,10 @@ void GSkyMap::read_wcs(const GFitsImage&amp; image)
2405 2580 alloc_pixels();
2406 2581  
2407 2582 // Read image
2408   - if (image.naxis() == 2) {
2409   - double* ptr = m_pixels;
2410   - for (int iy = 0; iy < m_num_y; ++iy) {
2411   - for (int ix = 0; ix < m_num_x; ++ix) {
2412   - *ptr++ = image.pixel(ix,iy);
2413   - }
2414   - }
2415   - }
2416   - else {
2417   - double* ptr = m_pixels;
2418   - for (int imap = 0; imap < m_num_maps; ++imap) {
2419   - for (int iy = 0; iy < m_num_y; ++iy) {
2420   - for (int ix = 0; ix < m_num_x; ++ix) {
2421   - *ptr++ = image.pixel(ix,iy,imap);
2422   - }
2423   - }
2424   - }
  2583 + double* ptr = m_pixels;
  2584 + int size = m_num_pixels * m_num_maps;
  2585 + for (int i = 0; i < size; ++i) {
  2586 + *ptr++ = (double)image.pixel(i);
2425 2587 }
2426 2588  
2427 2589 // Check if the map is too large
... ... @@ -2570,12 +2732,22 @@ GFitsImageDouble* GSkyMap::create_wcs_hdu(void) const
2570 2732 if (size > 0) {
2571 2733  
2572 2734 // Set axis parameters for image construction
2573   - int naxis = (m_num_maps == 1) ? 2 : 3;
2574   - int naxes[] = {m_num_x, m_num_y, m_num_maps};
  2735 + int naxis = (m_num_maps == 1) ? 2 : 2 + ndim();
  2736 +
  2737 + // Set dimensions of all axes
  2738 + int* naxes = new int[naxis];
  2739 + naxes[0] = m_num_x;
  2740 + naxes[1] = m_num_y;
  2741 + for (int i = 0; i < ndim(); ++i) {
  2742 + naxes[2+i] = m_shape[i];
  2743 + }
2575 2744  
2576 2745 // Allocate image
2577 2746 hdu = new GFitsImageDouble(naxis, naxes);
2578 2747  
  2748 + // Free dimensions of all axes
  2749 + delete [] naxes;
  2750 +
2579 2751 // Store data in image
2580 2752 if (naxis == 2) {
2581 2753 double* ptr = m_pixels;
... ...
test/test_GSky.cpp
  View file @ac6b92b
... ... @@ -842,6 +842,12 @@ void TestGSky::test_GSkyMap(void)
842 842 GSkyMap map_dst("CAR", "GAL", 0.0, 0.0, -0.1, 0.1, 100, 100, 2);
843 843 GSkyMap map_new("CAR", "GAL", 0.0, 0.0, -0.1, 0.1, 100, 100, 2);
844 844  
  845 + // Test map dimensions and shape
  846 + test_value(map_src.nmaps(), 2, "Check that sky map contains 2 maps");
  847 + test_value(map_src.ndim(), 1, "Check that sky map has one dimension");
  848 + test_value(map_src.shape().size(), 1, "Check that sky map has a shape size of 1");
  849 + test_value(map_src.shape()[0], 2, "Check that sky map has 2 maps");
  850 +
845 851 // Fill map pixels
846 852 double total_src = 0.0;
847 853 for (int pix = 0; pix < map_src.npix(); ++pix) {
... ... @@ -1027,6 +1033,50 @@ void TestGSky::test_GSkyMap(void)
1027 1033 test_value(total_extract, total_src, 1.0e-3, "Test extract() method with 2 maps");
1028 1034 test_value(map_extract.nmaps(), 2, "Test extract() method with 2 maps");
1029 1035  
  1036 + // Define one more map for shaping manipulation
  1037 + GSkyMap map_shape("CAR", "GAL", 0.0, 0.0, -1.0, 1.0, 10, 10, 12);
  1038 +
  1039 + // Test initial map dimensions and shape
  1040 + test_value(map_shape.nmaps(), 12, "Check that sky map contains 12 maps");
  1041 + test_value(map_shape.ndim(), 1, "Check that sky map has one dimension");
  1042 + test_value(map_shape.shape().size(), 1, "Check that sky map has a shape size of 1");
  1043 + test_value(map_shape.shape()[0], 12, "Check that sky map has 12 maps");
  1044 +
  1045 + // Set new map shape
  1046 + map_shape.shape(3,4);
  1047 +
  1048 + // Test map dimensions and shape
  1049 + test_value(map_shape.nmaps(), 12, "Check that sky map contains 12 maps");
  1050 + test_value(map_shape.ndim(), 2, "Check that sky map has two dimensions");
  1051 + test_value(map_shape.shape().size(), 2, "Check that sky map has a shape size of 2");
  1052 + test_value(map_shape.shape()[0], 3, "Check that sky map has 3 maps in first dimension");
  1053 + test_value(map_shape.shape()[1], 4, "Check that sky map has 4 maps in second dimension");
  1054 +
  1055 + // Set new map shape
  1056 + map_shape.shape(2,3,2);
  1057 +
  1058 + // Test map dimensions and shape
  1059 + test_value(map_shape.nmaps(), 12, "Check that sky map contains 12 maps");
  1060 + test_value(map_shape.ndim(), 3, "Check that sky map has three dimensions");
  1061 + test_value(map_shape.shape().size(), 3, "Check that sky map has a shape size of 3");
  1062 + test_value(map_shape.shape()[0], 2, "Check that sky map has 2 maps in first dimension");
  1063 + test_value(map_shape.shape()[1], 3, "Check that sky map has 3 maps in second dimension");
  1064 + test_value(map_shape.shape()[2], 2, "Check that sky map has 2 maps in third dimension");
  1065 +
  1066 + // Save map
  1067 + map_shape.save("test_map_shape.fits", true);
  1068 +
  1069 + // Load map
  1070 + GSkyMap map_load_shape("test_map_shape.fits");
  1071 +
  1072 + // Test map dimensions and shape
  1073 + test_value(map_load_shape.nmaps(), 12, "Check that sky map contains 12 maps");
  1074 + test_value(map_load_shape.ndim(), 3, "Check that sky map has three dimensions");
  1075 + test_value(map_load_shape.shape().size(), 3, "Check that sky map has a shape size of 3");
  1076 + test_value(map_load_shape.shape()[0], 2, "Check that sky map has 2 maps in first dimension");
  1077 + test_value(map_load_shape.shape()[1], 3, "Check that sky map has 3 maps in second dimension");
  1078 + test_value(map_load_shape.shape()[2], 2, "Check that sky map has 2 maps in third dimension");
  1079 +
1030 1080 // Exit test
1031 1081 return;
1032 1082 }
... ...
test/test_GSky.py
  View file @ac6b92b
1 1 # ==========================================================================
2 2 # This module performs unit tests for the GammaLib sky module.
3 3 #
4   -# Copyright (C) 2012-2015 Juergen Knoedlseder
  4 +# Copyright (C) 2012-2016 Juergen Knoedlseder
5 5 #
6 6 # This program is free software: you can redistribute it and/or modify
7 7 # it under the terms of the GNU General Public License as published by
... ... @@ -71,6 +71,8 @@ class Test(gammalib.GPythonTestSuite):
71 71 self.name("sky")
72 72  
73 73 # Append tests
  74 + self.append(self.test_methods, "Test sky map methods")
  75 + self.append(self.test_operators, "Test sky map operators")
74 76 self.append(self.test_skymap_healpix, "Test HEALPix map")
75 77 self.append(self.test_skymap_ait, "Test AIT projection map")
76 78 self.append(self.test_skymap_azp, "Test AZP projection map")
... ... @@ -80,11 +82,109 @@ class Test(gammalib.GPythonTestSuite):
80 82 self.append(self.test_skymap_stg, "Test STG projection map")
81 83 self.append(self.test_skymap_tan, "Test TAN projection map")
82 84 self.append(self.test_fk5_to_galactic, "Test FK5 to Galactic coordinate conversion")
83   - self.append(self.test_operators, "Test skymap operators")
84 85  
85 86 # Return
86 87 return
87 88  
  89 + # Test sky map methods
  90 + def test_methods(self):
  91 + """
  92 + Test the sky map methods.
  93 + """
  94 + # Setup sky maps
  95 + map = gammalib.GSkyMap("CAR", "CEL", 83.63, 22.01, -3.7, 2.6, 10, 8, 12)
  96 +
  97 + # Set sky map shape
  98 + map.shape(3,4)
  99 +
  100 + # Test map dimensions and shape
  101 + self.test_value(map.npix(), 80, "Check that sky map has 80 pixels");
  102 + self.test_value(map.nx(), 10, "Check that sky map has X dimension of 10");
  103 + self.test_value(map.ny(), 8, "Check that sky map has Y dimension of 8");
  104 + self.test_value(map.nmaps(), 12, "Check that sky map contains 12 maps");
  105 + self.test_value(map.ndim(), 2, "Check that sky map has two dimensions");
  106 + self.test_value(len(map.shape()), 2, "Check that sky map has a shape size of 2");
  107 + self.test_value(map.shape()[0], 3, "Check that sky map has 3 maps in first dimension");
  108 + self.test_value(map.shape()[1], 4, "Check that sky map has 4 maps in second dimension");
  109 +
  110 + # Set new map shape
  111 + map.shape(2,3,2);
  112 +
  113 + # Test map dimensions and shape
  114 + self.test_value(map.nmaps(), 12, "Check that sky map contains 12 maps");
  115 + self.test_value(map.ndim(), 3, "Check that sky map has three dimensions");
  116 + self.test_value(len(map.shape()), 3, "Check that sky map has a shape size of 3");
  117 + self.test_value(map.shape()[0], 2, "Check that sky map has 2 maps in first dimension");
  118 + self.test_value(map.shape()[1], 3, "Check that sky map has 3 maps in second dimension");
  119 + self.test_value(map.shape()[2], 2, "Check that sky map has 2 maps in third dimension");
  120 +
  121 + # Return
  122 + return
  123 +
  124 + # Test sky map operators
  125 + def test_operators(self):
  126 + """
  127 + Test the sky map operators.
  128 + """
  129 + # Setup sky maps
  130 + map = gammalib.GSkyMap("CAR", "CEL", 83.6331, 22.0145, -3.7, 2.6, 2, 2)
  131 + map[0] = 1.0
  132 + map[1] = 2.0
  133 + map[2] = 3.0
  134 + map[3] = 4.0
  135 + map_b = map.copy()
  136 +
  137 + # Addition operator
  138 + map += map_b
  139 + self.test_value(map[0], 2.0)
  140 + self.test_value(map[1], 4.0)
  141 + self.test_value(map[2], 6.0)
  142 + self.test_value(map[3], 8.0)
  143 +
  144 + # Multiplication operator
  145 + map *= map_b
  146 + self.test_value(map[0], 2.0)
  147 + self.test_value(map[1], 8.0)
  148 + self.test_value(map[2], 18.0)
  149 + self.test_value(map[3], 32.0)
  150 +
  151 + # Subtraction operator
  152 + map -= map_b
  153 + self.test_value(map[0], 1.0)
  154 + self.test_value(map[1], 6.0)
  155 + self.test_value(map[2], 15.0)
  156 + self.test_value(map[3], 28.0)
  157 +
  158 + # Division operator
  159 + map /= map_b
  160 + self.test_value(map[0], 1.0)
  161 + self.test_value(map[1], 3.0)
  162 + self.test_value(map[2], 5.0)
  163 + self.test_value(map[3], 7.0)
  164 +
  165 + # Scaling operator
  166 + map *= 2.0
  167 + self.test_value(map[0], 2.0)
  168 + self.test_value(map[1], 6.0)
  169 + self.test_value(map[2], 10.0)
  170 + self.test_value(map[3], 14.0)
  171 +
  172 + # Division operator
  173 + map /= 2.0
  174 + self.test_value(map[0], 1.0)
  175 + self.test_value(map[1], 3.0)
  176 + self.test_value(map[2], 5.0)
  177 + self.test_value(map[3], 7.0)
  178 +
  179 + # Access operator (tests also proper iteration)
  180 + sum = 0.0
  181 + for pix in map:
  182 + sum += pix
  183 + self.test_value(sum, 16.0)
  184 +
  185 + # Return
  186 + return
  187 +
88 188 # Generic skymap pixel transformation test
89 189 def test_skymap_pixels(self, pixels, map):
90 190 """
... ... @@ -315,68 +415,3 @@ class Test(gammalib.GPythonTestSuite):
315 415  
316 416 # Return
317 417 return
318   -
319   - # Test skymap operators
320   - def test_operators(self):
321   - """
322   - Test the skymap operators.
323   - """
324   - # Setup skymaps
325   - map = gammalib.GSkyMap("CAR", "CEL", 83.6331, 22.0145, -3.7, 2.6, 2, 2)
326   - map[0] = 1.0
327   - map[1] = 2.0
328   - map[2] = 3.0
329   - map[3] = 4.0
330   - map_b = map.copy()
331   -
332   - # Addition operator
333   - map += map_b
334   - self.test_value(map[0], 2.0)
335   - self.test_value(map[1], 4.0)
336   - self.test_value(map[2], 6.0)
337   - self.test_value(map[3], 8.0)
338   -
339   - # Multiplication operator
340   - map *= map_b
341   - self.test_value(map[0], 2.0)
342   - self.test_value(map[1], 8.0)
343   - self.test_value(map[2], 18.0)
344   - self.test_value(map[3], 32.0)
345   -
346   - # Subtraction operator
347   - map -= map_b
348   - self.test_value(map[0], 1.0)
349   - self.test_value(map[1], 6.0)
350   - self.test_value(map[2], 15.0)
351   - self.test_value(map[3], 28.0)
352   -
353   - # Division operator
354   - map /= map_b
355   - self.test_value(map[0], 1.0)
356   - self.test_value(map[1], 3.0)
357   - self.test_value(map[2], 5.0)
358   - self.test_value(map[3], 7.0)
359   -
360   - # Scaling operator
361   - map *= 2.0
362   - self.test_value(map[0], 2.0)
363   - self.test_value(map[1], 6.0)
364   - self.test_value(map[2], 10.0)
365   - self.test_value(map[3], 14.0)
366   -
367   - # Division operator
368   - map /= 2.0
369   - self.test_value(map[0], 1.0)
370   - self.test_value(map[1], 3.0)
371   - self.test_value(map[2], 5.0)
372   - self.test_value(map[3], 7.0)
373   -
374   - # Access operator (tests also proper iteration)
375   - sum = 0.0
376   - for pix in map:
377   - sum += pix
378   - self.test_value(sum, 16.0)
379   -
380   - # Return
381   - return
382   -
383 418 \ No newline at end of file
... ...