GammaLib  2.1.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GSkyMap.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GSkyMap.cpp - Sky map class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2010-2023 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GSkyMap.cpp
23  * @brief Sky map class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GException.hpp"
32 #include "GMath.hpp"
33 #include "GTools.hpp"
34 #include "GFft.hpp"
35 #include "GFits.hpp"
36 #include "GFitsTable.hpp"
37 #include "GFitsBinTable.hpp"
38 #include "GFitsImage.hpp"
39 #include "GFitsImageDouble.hpp"
40 #include "GFitsTableDoubleCol.hpp"
41 #include "GSkyDirs.hpp"
42 #include "GSkyMap.hpp"
43 #include "GHealpix.hpp"
44 #include "GWcsRegistry.hpp"
45 #include "GWcs.hpp"
46 #include "GMatrix.hpp"
47 #include "GVector.hpp"
48 #include "GVOClient.hpp"
49 #include "GSkyRegion.hpp"
50 #include "GSkyRegionCircle.hpp"
51 #include "GSkyRegions.hpp"
52 
53 /* __ Method name definitions ____________________________________________ */
54 #define G_CONSTRUCT_HPX "GSkyMap::GSkyMap(std::string&, int&,"\
55  " std::string&, int&)"
56 #define G_CONSTRUCT_MAP "GSkyMap::GSkyMap(std::string&, std::string&,"\
57  " double&, double&, double& double&, int&, int&, int&)"
58 #define G_NMAPS "GSkyMap::nmaps(int&)"
59 #define G_SHAPE "GSkyMap::shape(std::vector<int>&)"
60 #define G_OP_UNARY_ADD "GSkyMap::operator+=(GSkyMap&)"
61 #define G_OP_UNARY_SUB "GSkyMap::operator-=(GSkyMap&)"
62 #define G_OP_UNARY_MUL "GSkyMap::operator-=(GSkyMap&)"
63 #define G_OP_UNARY_DIV "GSkyMap::operator/=(GSkyMap&)"
64 #define G_OP_UNARY_DIV2 "GSkyMap::operator/=(double&)"
65 #define G_OP_ACCESS_1D "GSkyMap::operator(int&, int&)"
66 #define G_OP_ACCESS_2D "GSkyMap::operator(GSkyPixel&, int&)"
67 #define G_OP_VALUE "GSkyMap::operator(GSkyDir&, int&)"
68 #define G_INX2DIR "GSkyMap::inx2dir(int&)"
69 #define G_PIX2DIR "GSkyMap::pix2dir(GSkyPixel&)"
70 #define G_DIR2INX "GSkyMap::dir2inx(GSkyDir&)"
71 #define G_DIR2PIX "GSkyMap::dir2pix(GSkyDir&)"
72 #define G_FLUX1 "GSkyMap::flux(int&)"
73 #define G_FLUX2 "GSkyMap::flux(GSkyPixel&)"
74 #define G_FLUX3 "GSkyMap::flux(GSkyRegion&, int&)"
75 #define G_FLUX4 "GSkyMap::flux(GSkyRegions&, int&)"
76 #define G_SOLIDANGLE1 "GSkyMap::solidangle(int&)"
77 #define G_SOLIDANGLE2 "GSkyMap::solidangle(GSkyPixel&)"
78 #define G_OVERLAPS "GSkyMap::overlaps(GSkyRegion&)"
79 #define G_EXTRACT "GSkyMap::extract(int&, int&)"
80 #define G_EXTRACT_INT "GSkyMap::extract(int&, int&, int&, int&)"
81 #define G_EXTRACT_REG "GSkyMap::extract(GSkyRegions&)"
82 #define G_READ "GSkyMap::read(GFitsHDU&)"
83 #define G_SET_WCS "GSkyMap::set_wcs(std::string&, std::string&, double&,"\
84  " double&, double&, double&, double&, double&)"
85 #define G_READ_HEALPIX "GSkyMap::read_healpix(GFitsTable*)"
86 #define G_READ_WCS "GSkyMap::read_wcs(GFitsImage*)"
87 #define G_ALLOC_WCS "GSkyMap::alloc_wcs(GFitsImage*)"
88 #define G_SQRT "GSkyMap sqrt(GSkyMap&)"
89 #define G_LOG "GSkyMap log(GSkyMap&)"
90 #define G_LOG10 "GSkyMap log10(GSkyMap&)"
91 #define G_CONVOLUTION_KERNEL "GSkyMap::convolution_kernel(std::string&, "\
92  "double&, bool&)"
93 
94 /* __ Macros _____________________________________________________________ */
95 
96 /* __ Coding definitions _________________________________________________ */
97 
98 /* __ Debug definitions __________________________________________________ */
99 //#define G_READ_HEALPIX_DEBUG // Debug read_healpix
100 //#define G_READ_WCS_DEBUG // Debug read_wcs
101 
102 /* __ Prototype __________________________________________________________ */
103 
104 
105 /*==========================================================================
106  = =
107  = Constructors/destructors =
108  = =
109  ==========================================================================*/
110 
111 /***********************************************************************//**
112  * @brief Void constructor
113  *
114  * Constructs an empty sky map.
115  ***************************************************************************/
117 {
118  // Initialise class members for clean destruction
119  init_members();
120 
121  // Return
122  return;
123 }
124 
125 
126 /***********************************************************************//**
127  * @brief FITS file constructor
128  *
129  * @param[in] filename FITS file name.
130  *
131  * Constructs a sky map by loading data from a FITS file. See the load()
132  * method for more information about the FITS formats that are supported.
133  ***************************************************************************/
134 GSkyMap::GSkyMap(const GFilename& filename)
135 {
136  // Initialise class members for clean destruction
137  init_members();
138 
139  // Load skymap
140  load(filename);
141 
142  // Return
143  return;
144 }
145 
146 
147 /***********************************************************************//**
148  * @brief FITS HDU constructor
149  *
150  * @param[in] hdu FITS HDU.
151  *
152  * Constructs a sky map by fetching data from a FITS HDU. See the read()
153  * method for more information.
154  ***************************************************************************/
156 {
157  // Initialise class members for clean destruction
158  init_members();
159 
160  // Read skymap from HDU
161  read(hdu);
162 
163  // Return
164  return;
165 }
166 
167 
168 /***********************************************************************//**
169  * @brief Healpix sky map constructor
170  *
171  * @param[in] coords Coordinate System (CEL or GAL).
172  * @param[in] nside Nside parameter.
173  * @param[in] order Pixel ordering (RING or NEST).
174  * @param[in] nmaps Number of maps in set.
175  *
176  * @exception GException::invalid_argument
177  * Invalid @p nmaps parameter.
178  *
179  * Constructs @p nmaps identical all sky maps in Healpix pixelisation. All
180  * pixels of the sky maps will be initialised to values of zero.
181  ***************************************************************************/
182 GSkyMap::GSkyMap(const std::string& coords,
183  const int& nside,
184  const std::string& order,
185  const int& nmaps)
186 {
187  // Initialise class members for clean destruction
188  init_members();
189 
190  // Check if nmaps parameter is >0
191  if (nmaps < 1) {
192  std::string msg = "Number of maps "+gammalib::str(nmaps)+" is smaller "
193  "than 1. Please specify at least one map.";
195  }
196 
197  // Allocate Healpix projection
198  GHealpix* projection = new GHealpix(nside, order, coords);
199  m_proj = projection;
200 
201  // Set number of pixels, number of maps and shape of maps
202  m_num_pixels = projection->npix();
203  m_num_maps = nmaps;
204  m_shape.push_back(nmaps);
205 
206  // Allocate pixels
208 
209  // Return
210  return;
211 }
212 
213 
214 /***********************************************************************//**
215  * @brief WCS sky map constructor
216  *
217  * @param[in] wcs World Coordinate System.
218  * @param[in] coords Coordinate System (CEL or GAL).
219  * @param[in] x X coordinate of sky map centre (deg).
220  * @param[in] y Y coordinate of sky map centre (deg).
221  * @param[in] dx Pixel size in x direction at centre (deg/pixel).
222  * @param[in] dy Pixel size in y direction at centre (deg/pixel).
223  * @param[in] nx Number of pixels in x direction.
224  * @param[in] ny Number of pixels in y direction.
225  * @param[in] nmaps Number of maps in set (default=1).
226  *
227  * @exception GException::invalid_argument
228  * Invalid number of pixels or maps.
229  *
230  * Constructs @p nmaps identical all sky maps in World Coordinate System
231  * projection. All pixels of the sky maps will be initialised to values of
232  * zero.
233  ***************************************************************************/
234 GSkyMap::GSkyMap(const std::string& wcs,
235  const std::string& coords,
236  const double& x,
237  const double& y,
238  const double& dx,
239  const double& dy,
240  const int& nx,
241  const int& ny,
242  const int& nmaps)
243 {
244  // Initialise class members for clean destruction
245  init_members();
246 
247  // Check parameters
248  if (nx < 1) {
249  std::string msg = "Number of pixels "+gammalib::str(nx)+" in x "
250  "direction is smaller than 1. Please specify at "
251  "least one pixel in x direction.";
253  }
254  if (ny < 1) {
255  std::string msg = "Number of pixels "+gammalib::str(ny)+" in y "
256  "direction is smaller than 1. Please specify at "
257  "least one pixel in y direction.";
259  }
260  if (nmaps < 1) {
261  std::string msg = "Number of maps "+gammalib::str(nmaps)+" is smaller "
262  "than 1. Please specify at least one map.";
264  }
265 
266  // Set WCS
267  double crval1 = x;
268  double crval2 = y;
269  double crpix1 = double(nx+1)/2.0;
270  double crpix2 = double(ny+1)/2.0;
271  double cdelt1 = dx;
272  double cdelt2 = dy;
273  set_wcs(wcs, coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2);
274 
275  // Set number of pixels, number of maps and shape of maps
276  m_num_x = nx;
277  m_num_y = ny;
279  m_num_maps = nmaps;
280  m_shape.push_back(nmaps);
281 
282  // Allocate pixels
283  m_pixels = GNdarray(nx, ny, nmaps);
284 
285  // Return
286  return;
287 }
288 
289 
290 /***********************************************************************//**
291  * @brief Copy constructor
292  *
293  * @param[in] map Sky map.
294  *
295  * Constructs sky maps by copying data from another sky map object.
296  ***************************************************************************/
298 {
299  // Initialise class members for clean destruction
300  init_members();
301 
302  // Copy members
303  copy_members(map);
304 
305  // Return
306  return;
307 }
308 
309 
310 /***********************************************************************//**
311  * @brief Destructor
312  ***************************************************************************/
314 {
315  // Free members
316  free_members();
317 
318  // Return
319  return;
320 }
321 
322 
323 /*==========================================================================
324  = =
325  = Operators =
326  = =
327  ==========================================================================*/
328 
329 /***********************************************************************//**
330  * @brief Assignment operator
331  *
332  * @param[in] map Sky map.
333  * @return Sky map.
334  *
335  * Assigns one sky map to another.
336  ***************************************************************************/
338 {
339  // Execute only if object is not identical
340  if (this != &map) {
341 
342  // Free members
343  free_members();
344 
345  // Initialise private members for clean destruction
346  init_members();
347 
348  // Copy members
349  copy_members(map);
350 
351  } // endif: object was not identical
352 
353  // Return this object
354  return *this;
355 }
356 
357 
358 /***********************************************************************//**
359  * @brief Value setting operator
360  *
361  * @param[in] value Value.
362  * @return Sky map.
363  *
364  * Sets all sky map pixels to the specified @p value.
365  ***************************************************************************/
366 GSkyMap& GSkyMap::operator=(const double& value)
367 {
368  // Get number of pixels
369  int num = m_pixels.size();
370 
371  // Loop over all pixels
372  for (int i = 0; i < num; ++i) {
373  m_pixels(i) = value;
374  }
375 
376  // Return this object
377  return *this;
378 }
379 
380 
381 /***********************************************************************//**
382  * @brief Map addition operator
383  *
384  * @param[in] map Sky map.
385  * @return Sky map.
386  *
387  * @exception GException::invalid_value
388  * Mismatch between number of maps in skymap object.
389  *
390  * Adds the content of @p map to the skymap. The operator only works on sky
391  * maps with an identical number of maps.
392  *
393  * If the sky maps have the same definitions, the operator will add the
394  * values of the sky map pixels. Otherwise, the map values are added by
395  * bi-linearily interpolating the values in the source sky map, allowing thus
396  * for a reprojection of sky map values.
397  ***************************************************************************/
399 {
400  // Check if number of layers are identical
401  if (map.nmaps() != nmaps()) {
402  std::string msg = "Mismatch of number of maps in skymap object"
403  " ("+gammalib::str(nmaps())+" maps in destination"
404  " map, "+gammalib::str(map.nmaps())+" in source"
405  " map).";
407  }
408 
409  // If maps are identical then add the pixels
410  if (is_same(map)) {
411 
412  // Loop over all pixels of sky map
413  for (int index = 0; index < npix(); ++index) {
414 
415  // Loop over all layers
416  for (int layer = 0; layer < nmaps(); ++layer) {
417 
418  // Add pixel
419  (*this)(index, layer) += map(index, layer);
420 
421  } // endfor: looped over all layers
422 
423  } // endfor: looped over all pixels
424 
425  }
426 
427  // ... otherwise use interpolation method and add map
428  else {
429 
430  // Loop over all pixels of sky map
431  for (int index = 0; index < npix(); ++index) {
432 
433  // Get sky direction of actual pixel
434  GSkyDir dir = inx2dir(index);
435 
436  // Loop over all layers
437  for (int layer = 0; layer < nmaps(); ++layer) {
438 
439  // Add value
440  (*this)(index, layer) += map(dir, layer);
441 
442  } // endfor: looped over all layers
443 
444  } // endfor: looped over all pixels
445 
446  } // endelse: interpolation method
447 
448  // Return this object
449  return *this;
450 }
451 
452 
453 /***********************************************************************//**
454  * @brief Value addition operator
455  *
456  * @param[in] value Value.
457  * @return Sky map.
458  *
459  * Add @p value to all sky map pixels.
460  ***************************************************************************/
461 GSkyMap& GSkyMap::operator+=(const double& value)
462 {
463  // Add value to all pixels
464  m_pixels += value;
465 
466  // Return this object
467  return *this;
468 }
469 
470 
471 /***********************************************************************//**
472  * @brief Map subtraction operator
473  *
474  * @param[in] map Sky map.
475  * @return Sky map.
476  *
477  * @exception GException::invalid_value
478  * Mismatch between number of maps in skymap object.
479  *
480  * Subtracts the content of @p map from the skymap. The operator only works
481  * on sky maps with an identical number of layers.
482  *
483  * If the sky maps have the same definitions, the operator will subtract the
484  * values of the sky map pixels. Otherwise, the map values are subtracted by
485  * bi-linearily interpolating the values in the source sky map, allowing thus
486  * for a reprojection of sky map values.
487  ***************************************************************************/
489 {
490  // Check if number of layers are identical
491  if (map.nmaps() != nmaps()) {
492  std::string msg = "Mismatch of number of maps in skymap object"
493  " ("+gammalib::str(nmaps())+" maps in destination"
494  " map, "+gammalib::str(map.nmaps())+" in source"
495  " map).";
497  }
498 
499  // If maps are identical then subtract the pixels
500  if (is_same(map)) {
501 
502  // Loop over all pixels of sky map
503  for (int index = 0; index < npix(); ++index) {
504 
505  // Loop over all layers
506  for (int layer = 0; layer < nmaps(); ++layer) {
507 
508  // Subtract pixel
509  (*this)(index, layer) -= map(index, layer);
510 
511  } // endfor: looped over all layers
512 
513  } // endfor: looped over all pixels
514 
515  }
516 
517  // ... otherwise use interpolation method and subtract map
518  else {
519 
520  // Loop over all pixels of sky map
521  for (int index = 0; index < npix(); ++index) {
522 
523  // Get sky direction of actual pixel
524  GSkyDir dir = inx2dir(index);
525 
526  // Loop over all layers
527  for (int layer = 0; layer < nmaps(); ++layer) {
528 
529  // Subtract value
530  (*this)(index, layer) -= map(dir, layer);
531 
532  } // endfor: looped over all layers
533 
534  } // endfor: looped over all pixels
535 
536  } // endelse: interpolation method
537 
538  // Return this object
539  return *this;
540 }
541 
542 
543 /***********************************************************************//**
544  * @brief Value subtraction operator
545  *
546  * @param[in] value Value.
547  * @return Sky map.
548  *
549  * Subtracts @p value from all sky map pixels.
550  ***************************************************************************/
551 GSkyMap& GSkyMap::operator-=(const double& value)
552 {
553  // Subtract value from all pixels
554  m_pixels -= value;
555 
556  // Return this object
557  return *this;
558 }
559 
560 
561 /***********************************************************************//**
562  * @brief Multiplication operator
563  *
564  * @param[in] map Sky map.
565  * @return Sky map.
566  *
567  * @exception GException::invalid_value
568  * Mismatch between number of maps in skymap object.
569  *
570  * Multiplies the content of @p map from the skymap. The operator only works
571  * on sky maps with an identical number of layers.
572  *
573  * If the sky maps have the same definitions, the operator will multiply the
574  * values of the sky map pixels. Otherwise, the map values are multiplied by
575  * bi-linearily interpolating the values in the source sky map, allowing thus
576  * for a reprojection of sky map values.
577  ***************************************************************************/
579 {
580  // Check if number of layers are identical
581  if (map.nmaps() != nmaps()) {
582  std::string msg = "Mismatch of number of maps in skymap object"
583  " ("+gammalib::str(nmaps())+" maps in destination"
584  " map, "+gammalib::str(map.nmaps())+" in source"
585  " map).";
587  }
588 
589  // If maps are identical then multiply the pixels
590  if (is_same(map)) {
591 
592  // Loop over all pixels of sky map
593  for (int index = 0; index < npix(); ++index) {
594 
595  // Loop over all layers
596  for (int layer = 0; layer < nmaps(); ++layer) {
597 
598  // Multiply pixel
599  (*this)(index, layer) *= map(index, layer);
600 
601  } // endfor: looped over all layers
602 
603  } // endfor: looped over all pixels
604 
605  }
606 
607  // ... otherwise use interpolation method and multiply map
608  else {
609 
610  // Loop over all pixels of sky map
611  for (int index = 0; index < npix(); ++index) {
612 
613  // Get sky direction of actual pixel
614  GSkyDir dir = inx2dir(index);
615 
616  // Loop over all layers
617  for (int layer = 0; layer < nmaps(); ++layer) {
618 
619  // Multiply value
620  (*this)(index, layer) *= map(dir, layer);
621 
622  } // endfor: looped over all layers
623 
624  } // endfor: looped over all pixels
625 
626  } // endelse: interpolation method
627 
628  // Return this object
629  return *this;
630 }
631 
632 
633 /***********************************************************************//**
634  * @brief Sky map scaling operator
635  *
636  * @param[in] factor Scale factor.
637  * @return Sky map.
638  *
639  * Multiplies all pixels of the sky map by the given scale @p factor.
640  ***************************************************************************/
641 GSkyMap& GSkyMap::operator*=(const double& factor)
642 {
643  // Scale all pixels
644  m_pixels *= factor;
645 
646  // Return this object
647  return *this;
648 }
649 
650 
651 /***********************************************************************//**
652  * @brief Division operator
653  *
654  * @param[in] map Sky map.
655  * @return Sky map.
656  *
657  * @exception GException::invalid_value
658  * Mismatch between number of maps in skymap object.
659  *
660  * Divides the content of the actual skymap by the skymap @p map. The operator
661  * only works on sky maps with an identical number of layers.
662  *
663  * If the sky maps have the same definitions, the operator will divide the
664  * values of the sky map pixels. Otherwise, the map values are divided by
665  * bi-linearily interpolating the values in the source sky map, allowing thus
666  * for a reprojection of sky map values.
667  *
668  * On return, all pixels in @p map that are zero are silently set to zero in
669  * the skymap, hence no division by zero exception will be raised in case
670  * that zero values are encountered.
671  ***************************************************************************/
673 {
674  // Check if number of layers are identical
675  if (map.nmaps() != nmaps()) {
676  std::string msg = "Mismatch of number of maps in skymap object"
677  " ("+gammalib::str(nmaps())+" maps in destination"
678  " map, "+gammalib::str(map.nmaps())+" in source"
679  " map).";
681  }
682 
683  // If maps are identical then divide the pixels
684  if (is_same(map)) {
685 
686  // Loop over all pixels of sky map
687  for (int index = 0; index < npix(); ++index) {
688 
689  // Loop over all layers
690  for (int layer = 0; layer < nmaps(); ++layer) {
691 
692  // Get pixel value of the map by which the division will be
693  // done.
694  double value = map(index, layer);
695 
696  // Divide destination pixel by value. In case of a division
697  // by zero set the destination pixel also to zero.
698  if (value == 0.0) {
699  (*this)(index, layer) = 0.0;
700  }
701  else {
702  (*this)(index, layer) /= value;
703  }
704 
705  } // endfor: looped over all layers
706 
707  } // endfor: looped over all pixels
708 
709  }
710 
711  // ... otherwise use interpolation method and divide map
712  else {
713 
714  // Loop over all pixels of destination sky map
715  for (int index = 0; index < npix(); ++index) {
716 
717  // Get sky direction of actual pixel
718  GSkyDir dir = inx2dir(index);
719 
720  // Loop over all layers
721  for (int layer = 0; layer < nmaps(); ++layer) {
722 
723  // Get map value of the map by which the division will be
724  // done. The map is accessed using the sky direction, hence
725  // interpolating properly the map. If we're outside the map,
726  // zero is returned.
727  double value = map(dir, layer);
728 
729  // Divide destination pixel by value. In case of a division
730  // by zero set the destination pixel also to zero.
731  if (value == 0.0) {
732  (*this)(index, layer) = 0.0;
733  }
734  else {
735  (*this)(index, layer) /= value;
736  }
737 
738  } // endfor: looped over all layers
739 
740  } // endfor: looped over all pixels
741 
742  } // endelse: interpolation method
743 
744  // Return this object
745  return *this;
746 }
747 
748 
749 /***********************************************************************//**
750  * @brief Sky map division operator
751  *
752  * @param[in] factor Scale factor.
753  * @return Sky map.
754  *
755  * @exception GException::invalid_argument
756  * Division by zero error.
757  *
758  * Divides all pixels of the sky map by the given @p factor.
759  ***************************************************************************/
760 GSkyMap& GSkyMap::operator/=(const double& factor)
761 {
762  // Check for division by zero
763  if (factor == 0.0) {
764  std::string msg = "Trying to divide sky map pixels by zero.";
766  }
767 
768  // Divide all pixels by factor
769  m_pixels /= factor;
770 
771  // Return this object
772  return *this;
773 }
774 
775 
776 /***********************************************************************//**
777  * @brief Pixel index access operator
778  *
779  * @param[in] index Pixel index [0,...,npix()-1].
780  * @param[in] map Map index [0,...,nmaps()-1].
781  * @return Sky map pixel value.
782  *
783  * @exception GException::out_of_range
784  * Pixel index and/or map index are outside valid range.
785  *
786  * Access sky map pixel by its index, where the most quickly varying axis is
787  * the x axis of the map.
788  ***************************************************************************/
789 double& GSkyMap::operator()(const int& index, const int& map)
790 {
791  // Throw an error if pixel index or map index is not in valid range
792  #if defined(G_RANGE_CHECK)
793  if (index < 0 || index >= m_num_pixels) {
795  "Sky map pixel index",
796  index, m_num_pixels);
797  }
798  if (map < 0 || map >= m_num_maps) {
800  "Sky map map index",
801  map, m_num_maps);
802  }
803  #endif
804 
805  // Return reference to pixel value
806  return m_pixels(index+m_num_pixels*map);
807 }
808 
809 
810 /***********************************************************************//**
811  * @brief Pixel index access operator (const variant)
812  *
813  * @param[in] index Pixel index [0,...,npix()-1].
814  * @param[in] map Map index [0,...,nmaps()-1].
815  *
816  * @exception GException::out_of_range
817  * Pixel index and/or map index are outside valid range.
818  *
819  * Access sky map pixel by its index, where the most quickly varying axis is
820  * the x axis of the map.
821  ***************************************************************************/
822 const double& GSkyMap::operator()(const int& index, const int& map) const
823 {
824  // Throw an error if pixel index or map index is not in valid range
825  #if defined(G_RANGE_CHECK)
826  if (index < 0 || index >= m_num_pixels) {
828  "Sky map pixel index",
829  index, m_num_pixels);
830  }
831  if (map < 0 || map >= m_num_maps) {
833  "Sky map map index",
834  map, m_num_maps);
835  }
836  #endif
837 
838  // Return reference to pixel value
839  return m_pixels(index+m_num_pixels*map);
840 }
841 
842 
843 /***********************************************************************//**
844  * @brief Sky map pixel access operator
845  *
846  * @param[in] pixel Sky map pixel.
847  * @param[in] map Sky map index [0,...,nmaps()[.
848  *
849  * @exception GException::invalid_argument
850  * Sky pixel is not contained in sky map.
851  * @exception GException::out_of_range
852  * Sky map index is out of range.
853  *
854  * Access sky map pixel by its 2D index (x,y) that is implemented by the
855  * GSkyPixel class.
856  ***************************************************************************/
857 double& GSkyMap::operator()(const GSkyPixel& pixel, const int& map)
858 {
859  // Throw an error if pixel index or map index is not in valid range
860  #if defined(G_RANGE_CHECK)
861  if (!contains(pixel)) {
862  std::string msg = "Sky pixel ("+gammalib::str(pixel.x())+","+
863  gammalib::str(pixel.y())+") is not contained in "
864  "sky map comprised of ("+gammalib::str(m_num_x)+","+
865  gammalib::str(m_num_y)+") pixels. Please specify "
866  "a valid sky pixel.";
868  }
869  if (map < 0 || map >= m_num_maps) {
870  throw GException::out_of_range(G_OP_ACCESS_2D, "Sky map index",
871  map, m_num_maps);
872  }
873  #endif
874 
875  // Get pixel index
876  int index = pix2inx(pixel);
877 
878  // Return reference to pixel value
879  return m_pixels(index+m_num_pixels*map);
880 }
881 
882 
883 /***********************************************************************//**
884  * @brief Sky map pixel access operator
885  *
886  * @param[in] pixel Sky map pixel.
887  * @param[in] map Sky map index [0,...,nmaps()[.
888  *
889  * @exception GException::invalid_argument
890  * Sky pixel is not contained in sky map.
891  * @exception GException::out_of_range
892  * Sky map index is out of range.
893  *
894  * Access sky map pixel by its 2D index (x,y) that is implemented by the
895  * GSkyPixel class.
896  ***************************************************************************/
897 const double& GSkyMap::operator()(const GSkyPixel& pixel, const int& map) const
898 {
899  // Throw an error if pixel index or map index is not in valid range
900  #if defined(G_RANGE_CHECK)
901  if (!contains(pixel)) {
902  std::string msg = "Sky pixel ("+gammalib::str(pixel.x())+","+
903  gammalib::str(pixel.y())+") is not contained in "
904  "sky map comprised of ("+gammalib::str(m_num_x)+","+
905  gammalib::str(m_num_y)+") pixels. Please specify "
906  "a valid sky pixel.";
908  }
909  if (map < 0 || map >= m_num_maps) {
910  throw GException::out_of_range(G_OP_ACCESS_2D, "Sky map index",
911  map, m_num_maps);
912  }
913  #endif
914 
915  // Get pixel index
916  int index = pix2inx(pixel);
917 
918  // Return reference to pixel value
919  return m_pixels(index+m_num_pixels*map);
920 }
921 
922 
923 /***********************************************************************//**
924  * @brief Return interpolated skymap value for sky direction
925  *
926  * @param[in] dir Sky direction.
927  * @param[in] map Map index [0,...,nmaps()-1].
928  * @return Sky intensity.
929  *
930  * @exception GException::out_of_range
931  * Map index lies outside valid range.
932  *
933  * Returns the skymap value for a given sky direction, obtained by bi-linear
934  * interpolation of the neighbouring pixels. If the sky direction falls
935  * outside the area covered by the skymap, a value of 0 is returned.
936  *
937  * The method implements a computation cache that avoids recomputation of
938  * interpolation indices and weights in case that the same sky direction
939  * is requested several times. This speeds up use cases where skymap values
940  * for various map indices have to be returned for the same sky direction.
941  ***************************************************************************/
942 double GSkyMap::operator()(const GSkyDir& dir, const int& map) const
943 {
944  // Throw an error if the map index is not in valid range
945  #if defined(G_RANGE_CHECK)
946  if (map < 0 || map >= m_num_maps) {
948  "Sky map map index",
949  map, m_num_maps);
950  }
951  #endif
952 
953  // Initialise intensity
954  double intensity = 0.0;
955 
956  // If the pre-computation cache has not been set or if the sky direction
957  // has changed then compute the bilinear interpolator
958  if (!m_hascache || (dir != m_last_dir)) {
959 
960  // Perform computation for HealPix map
961  if (m_proj->size() == 1) {
962 
963  // Get Healpix interpolator
964  m_interpol = static_cast<GHealpix*>(m_proj)->interpolator(dir);
965 
966  // Set containment flag and store actual sky direction as last
967  // direction
968  m_hascache = true;
969  m_contained = true;
970  m_last_dir = dir;
971 
972  } // endif: we had a Healpix projection
973 
974  // ... otherwise perform computation for WCS map
975  else {
976 
977  // Determine sky pixel. At this point an exception may occur
978  // in case that the pixel cannot be represented by the
979  // relevant projection. We catch this exception here
980  try {
981 
982  // Determine sky pixel
983  GSkyPixel pixel = dir2pix(dir);
984 
985  // Set containment flag and store actual sky direction
986  // as last direction
987  m_hascache = true;
988  m_contained = contains(pixel);
989  m_last_dir = dir;
990 
991  // Continue only if pixel is within the map
992  if (m_contained) {
993 
994  // Signal that we have a wrap around in the x axis
995  double x_size = std::abs(m_num_x * static_cast<GWcs*>(m_proj)->cdelt(0));
996  bool x_wrap = (x_size > 359.99);
997 
998  // Get pixel index that is to the left-top of the actual
999  // pixel. We take care of the special case of negative pixel
1000  // indices which arise if we are in the first column or
1001  // first row of the map.
1002  int inx_x = int(pixel.x());
1003  int inx_y = int(pixel.y());
1004  if (pixel.x() < 0.0) {
1005  inx_x--;
1006  }
1007  if (pixel.y() < 0.0) {
1008  inx_y--;
1009  }
1010 
1011  // Set left and right indices for interpolation. The left
1012  // index needs to be non-negative and the right index needs
1013  // to be not larger that the number of pixels. We treat
1014  // here also the special case of wrap around in the x
1015  // axis that may occur if we have an allsky map.
1016  int inx_x_left = inx_x;
1017  int inx_y_left = inx_y;
1018  int inx_x_right = inx_x_left + 1;
1019  int inx_y_right = inx_y_left + 1;
1020  if (inx_x_left < 0) {
1021  if (x_wrap) {
1022  inx_x_left += m_num_x;
1023  }
1024  else {
1025  inx_x_left = 0;
1026  }
1027  }
1028  if (inx_x_right >= m_num_x) {
1029  if (x_wrap) {
1030  inx_x_right -= m_num_x;
1031  }
1032  else {
1033  inx_x_right = m_num_x - 1;
1034  }
1035  }
1036  if (inx_y_left < 0) {
1037  inx_y_left = 0;
1038  }
1039  if (inx_y_right >= m_num_y) {
1040  inx_y_right = m_num_y - 1;
1041  }
1042 
1043  // Set weighting factors for interpolation
1044  double wgt_x_right = (pixel.x() - inx_x);
1045  double wgt_x_left = 1.0 - wgt_x_right;
1046  double wgt_y_right = (pixel.y() - inx_y);
1047  double wgt_y_left = 1.0 - wgt_y_right;
1048 
1049  // Compute skymap pixel indices for bi-linear interpolation
1050  m_interpol.index1() = inx_x_left + inx_y_left * m_num_x;
1051  m_interpol.index2() = inx_x_left + inx_y_right * m_num_x;
1052  m_interpol.index3() = inx_x_right + inx_y_left * m_num_x;
1053  m_interpol.index4() = inx_x_right + inx_y_right * m_num_x;
1054 
1055  // Compute weighting factors for bi-linear interpolation
1056  m_interpol.weight1() = wgt_x_left * wgt_y_left;
1057  m_interpol.weight2() = wgt_x_left * wgt_y_right;
1058  m_interpol.weight3() = wgt_x_right * wgt_y_left;
1059  m_interpol.weight4() = wgt_x_right * wgt_y_right;
1060 
1061  } // endif: pixel was contained in map
1062 
1063  } // endtry: pixel computation was successful
1065  m_contained = false;
1066  }
1067 
1068  } // endelse: we had a WCS map
1069 
1070  } // endif: computation of the bi-linear interpolater was required
1071 
1072  // Compute the interpolated intensity if the pixel is contained in
1073  // the map
1074  if (m_contained) {
1075 
1076  // Compute map offset
1077  int offset = m_num_pixels * map;
1078 
1079  // Compute interpolated skymap value
1080  intensity = m_interpol(m_pixels.data()+offset);
1081 
1082  } // endif: direction was contained in map
1083 
1084  // Return intensity
1085  return intensity;
1086 }
1087 
1088 
1089 /*==========================================================================
1090  = =
1091  = Public methods =
1092  = =
1093  ==========================================================================*/
1094 
1095  /***********************************************************************//**
1096  * @brief Clear instance.
1097  *
1098  * Resets the sky map to the initial state.
1099  ***************************************************************************/
1100 void GSkyMap::clear(void)
1101 {
1102  // Free class members
1103  free_members();
1104 
1105  // Initialise members
1106  init_members();
1107 
1108  // Return
1109  return;
1110 }
1111 
1112 
1113 /***********************************************************************//**
1114  * @brief Clone sky map
1115  *
1116  * @return Pointer to deep copy of sky map.
1117  ***************************************************************************/
1119 {
1120  return new GSkyMap(*this);
1121 }
1122 
1123 
1124 /***********************************************************************//**
1125  * @brief Set number of maps
1126  *
1127  * @param[in] nmaps Number of maps.
1128  *
1129  * @exception GException::invalid_argument
1130  * Invalid number of maps specified.
1131  *
1132  * Redefines the number of maps in an GSkyMap object. If the number of maps
1133  * is increased with respect to the existing number, additional maps with
1134  * pixel values of zero are append to the object. Existing map pixel values
1135  * are kept. If the number of maps is decreased with respect to the existing
1136  * number, the excedent maps are dropped. The remaining map pixel values are
1137  * kept.
1138  ***************************************************************************/
1139 void GSkyMap::nmaps(const int& nmaps)
1140 {
1141  // Throw an exception if less than 1 map is required
1142  if (nmaps < 1) {
1143  std::string msg = "At least one map is required in an GSkyMap object."
1144  " Please specify an argument >=1.";
1146  }
1147 
1148  // If the map has pixels and if the number of maps changes then set a new
1149  // number of maps. Copy over any existing information.
1150  if (m_num_pixels > 0 && nmaps != m_num_maps) {
1151 
1152  // Allocate new array
1153  GNdarray new_pixels(m_num_pixels, nmaps);
1154 
1155  // Copy over existing pixels
1156  int num_copy = (nmaps > m_num_maps) ? m_num_maps : nmaps;
1157  num_copy *= m_num_pixels;
1158  for (int i = 0; i < num_copy; ++i) {
1159  new_pixels(i) = m_pixels(i);
1160  }
1161 
1162  // Set pixels to new array
1163  m_pixels = new_pixels;
1164 
1165  // Set number of maps
1166  m_num_maps = nmaps;
1167 
1168  // Set shape of maps to linear array since we do not know how to
1169  // arrange the new maps
1170  m_shape.clear();
1171  m_shape.push_back(nmaps);
1172 
1173  } // endif: map had pixels
1174 
1175  // Return
1176  return;
1177 }
1178 
1179 
1180 /***********************************************************************//**
1181  * @brief Set one-dimensional shape of maps
1182  *
1183  * @param[in] s1 Axis length of first dimension.
1184  ***************************************************************************/
1185 void GSkyMap::shape(const int& s1)
1186 {
1187  // Initialise vector
1188  std::vector<int> shape;
1189 
1190  // Set vector elements
1191  shape.push_back(s1);
1192 
1193  // Call vector method
1194  this->shape(shape);
1195 
1196  // Return
1197  return;
1198 }
1199 
1200 
1201 /***********************************************************************//**
1202  * @brief Set two-dimensional shape of maps
1203  *
1204  * @param[in] s1 Axis length of first dimension.
1205  * @param[in] s2 Axis length of second dimension.
1206  ***************************************************************************/
1207 void GSkyMap::shape(const int& s1, const int& s2)
1208 {
1209  // Initialise vector
1210  std::vector<int> shape;
1211 
1212  // Set vector elements
1213  shape.push_back(s1);
1214  shape.push_back(s2);
1215 
1216  // Call vector method
1217  this->shape(shape);
1218 
1219  // Return
1220  return;
1221 }
1222 
1223 
1224 /***********************************************************************//**
1225  * @brief Set three-dimensional shape of maps
1226  *
1227  * @param[in] s1 Axis length of first dimension.
1228  * @param[in] s2 Axis length of second dimension.
1229  * @param[in] s3 Axis length of second dimension.
1230  ***************************************************************************/
1231 void GSkyMap::shape(const int& s1, const int& s2, const int& s3)
1232 {
1233  // Initialise vector
1234  std::vector<int> shape;
1235 
1236  // Set vector elements
1237  shape.push_back(s1);
1238  shape.push_back(s2);
1239  shape.push_back(s3);
1240 
1241  // Call vector method
1242  this->shape(shape);
1243 
1244  // Return
1245  return;
1246 }
1247 
1248 
1249 /***********************************************************************//**
1250  * @brief Set shape of maps
1251  *
1252  * @param[in] shape Shape vector.
1253  *
1254  * @exception GException::invalid_argument
1255  * Invalid shape factorisation specified.
1256  *
1257  * Defines a shape for the maps in the object. The shape specifies how the
1258  * maps are arranged in a n-dimensional array.
1259  ***************************************************************************/
1260 void GSkyMap::shape(const std::vector<int>& shape)
1261 {
1262  // Computes the resulting number of maps
1263  int nmaps = 0;
1264  if (shape.size() > 0) {
1265  nmaps = 1;
1266  for (int i = 0; i < shape.size(); ++i) {
1267  nmaps *= shape[i];
1268  }
1269  }
1270 
1271  // Throw an exception if resulting number of maps is not equal to the
1272  // existing number of maps
1273  if (nmaps != m_num_maps) {
1274  std::string msg = "Sky map shaping requires "+gammalib::str(nmaps)+
1275  " maps while "+gammalib::str(m_num_maps)+" maps "
1276  "are available. Please specify an appropriate "
1277  "factorisation or modify the number of available "
1278  "maps.";
1280  }
1281 
1282  // Set shape
1283  m_shape = shape;
1284 
1285  // Return
1286  return;
1287 }
1288 
1289 
1290 /***********************************************************************//**
1291  * @brief Converts pixel index into sky map pixel
1292  *
1293  * @param[in] index Pixel index [0,...,npix()-1].
1294  * @return Sky map pixel.
1295  *
1296  * Converts the pixel @p index into a sky map pixel (GSkyPixel).The dimension
1297  * of GSkyPixel will be identical to the dimension of the sky map (i.e. a 1D
1298  * sky map leads to a 1D GSkyPixel object, a 2D sky map leads to a 2D
1299  * GSkyPixel object.
1300  ***************************************************************************/
1301 GSkyPixel GSkyMap::inx2pix(const int& index) const
1302 {
1303  // Initialise sky map pixel
1304  GSkyPixel pixel;
1305 
1306  // Get x and y indices
1307  if (m_num_x != 0) { //!< 2D sky map
1308  pixel.x(double(index % m_num_x));
1309  pixel.y(double(index / m_num_x));
1310  }
1311  else { //!< 1D sky map
1312  pixel.index(index);
1313  }
1314 
1315  // Return pixel
1316  return pixel;
1317 }
1318 
1319 
1320 /***********************************************************************//**
1321  * @brief Returns sky direction of pixel
1322  *
1323  * @param[in] index Pixel index [0,...,npix()-1].
1324  * @return Sky direction.
1325  *
1326  * @exception GException::invalid_value
1327  * No valid sky projection found.
1328  *
1329  * Returns sky direction for a given pixel index.
1330  ***************************************************************************/
1331 GSkyDir GSkyMap::inx2dir(const int& index) const
1332 {
1333  // Throw error if sky projection is not valid
1334  if (m_proj == NULL) {
1335  std::string msg = "Sky projection has not been defined.";
1337  }
1338 
1339  // Determine sky direction from pixel index.
1340  GSkyDir dir = m_proj->pix2dir(inx2pix(index));
1341 
1342  // Return sky direction
1343  return dir;
1344 }
1345 
1346 
1347 /***********************************************************************//**
1348  * @brief Returns sky direction of pixel
1349  *
1350  * @param[in] pixel Sky map pixel.
1351  * @return Sky direction.
1352  *
1353  * @exception GException::invalid_value
1354  * No valid sky projection found.
1355  * @exception GException::invalid_argument
1356  * 2D sky map pixel used to access 1D projection.
1357  *
1358  * Returns sky direction for a given sky map @p pixel.
1359  ***************************************************************************/
1361 {
1362  // Throw error if WCS is not valid
1363  if (m_proj == NULL) {
1364  std::string msg = "Sky projection has not been defined.";
1366  }
1367 
1368  // Initialise sky direction
1369  GSkyDir dir;
1370 
1371  // If pixel size matches the projection size then perform a straight
1372  // forward conversion
1373  if (m_proj->size() == pixel.size()) {
1374  dir = m_proj->pix2dir(pixel);
1375  }
1376 
1377  // ... otherwise, if we have a 2D projection but a 1D pixel then
1378  // interpret the pixel as the linear index in the pixel array
1379  else if (m_proj->size() == 2) {
1380  dir = m_proj->pix2dir(GSkyPixel(inx2pix(int(pixel))));
1381  }
1382 
1383  // ... otherwise we have a 1D projection but a 2D pixel. There is
1384  // no unambiguous way to handle this case, hence we throw an exception
1385  else {
1386  std::string msg = "A 2-dimensional sky map pixel "+pixel.print()+
1387  " is used to determine the sky direction for"
1388  " the 1-dimensional sky projection \""+
1389  m_proj->name()+"\"\n"
1390  "Please specify a 1-dimensional sky map pixel.";
1392  }
1393 
1394  // Return sky direction
1395  return dir;
1396 }
1397 
1398 
1399 /***********************************************************************//**
1400  * @brief Converts sky map pixel into pixel index
1401  *
1402  * @param[in] pixel Sky map pixel.
1403  * @return Pixel index [0,...,npix()-1].
1404  *
1405  * Converts a sky map @p pixel into the pixel index.
1406  ***************************************************************************/
1407 int GSkyMap::pix2inx(const GSkyPixel& pixel) const
1408 {
1409  // Initialise pixel index
1410  int index = 0;
1411 
1412  // Handle 1D sky map pixel
1413  if (pixel.is_1D()) {
1414  index = int(pixel);
1415  }
1416 
1417  // Handle 2D sky map pixel
1418  else if (pixel.is_2D()) {
1419 
1420  // Get x and y indices by rounding the (x,y) values
1421  int ix = int(pixel.x()+0.5);
1422  int iy = int(pixel.y()+0.5);
1423 
1424  // Set index
1425  index = ix + iy * m_num_x;
1426 
1427  }
1428 
1429  // Return index
1430  return index;
1431 }
1432 
1433 
1434 /***********************************************************************//**
1435  * @brief Returns pixel index for a given sky direction
1436  *
1437  * @param[in] dir Sky direction.
1438  * @return Pixel index [0,...,npix()-1].
1439  *
1440  * @exception GException::invalid_value
1441  * No valid sky projection found.
1442  *
1443  * Returns sky map pixel index for a given sky direction.
1444  ***************************************************************************/
1445 int GSkyMap::dir2inx(const GSkyDir& dir) const
1446 {
1447  // Throw error if WCS is not valid
1448  if (m_proj == NULL) {
1449  std::string msg = "Sky projection has not been defined.";
1451  }
1452 
1453  // Determine pixel index for a given sky direction
1454  int index = pix2inx(m_proj->dir2pix(dir));
1455 
1456  // Return pixel index
1457  return index;
1458 }
1459 
1460 
1461 /***********************************************************************//**
1462  * @brief Returns sky map pixel for a given sky direction
1463  *
1464  * @param[in] dir Sky direction.
1465  * @return Sky map pixel.
1466  *
1467  * @exception GException::invalid_value
1468  * No valid sky projection found.
1469  *
1470  * Returns sky map pixel for a given sky direction.
1471  ***************************************************************************/
1473 {
1474  // Throw error if WCS is not valid
1475  if (m_proj == NULL) {
1476  std::string msg = "Sky projection has not been defined.";
1478  }
1479 
1480  // Determine pixel for a given sky direction
1481  GSkyPixel pixel = m_proj->dir2pix(dir);
1482 
1483  // Return pixel
1484  return pixel;
1485 }
1486 
1487 
1488 /***********************************************************************//**
1489  * @brief Returns array with total number of counts for count maps
1490  *
1491  * @return Array of pixel counts sums.
1492  *
1493  * For a set of n count maps, the method returns a 1D array with n entries
1494  * that correspond to the sum of the pixel counts in each map.
1495  ***************************************************************************/
1497 {
1498  // Initialise output array
1499  GNdarray sum_array = GNdarray(m_num_maps);
1500 
1501  // Loop over maps and store sums in output array
1502  for (int i = 0; i < m_num_maps; ++i) {
1503 
1504  // Initialise sum for this map
1505  double sum = 0.0;
1506 
1507  // Compute map sum
1508  for (int k = 0; k < m_num_pixels; ++k) {
1509  sum += (*this)(k,i);
1510  }
1511 
1512  // Store map sum
1513  sum_array(i) = sum;
1514 
1515  } // endfor: looped over all maps
1516 
1517  // Return sum array
1518  return sum_array;
1519 }
1520 
1521 
1522 /***********************************************************************//**
1523  * @brief Returns flux in pixel
1524  *
1525  * @param[in] index Pixel index [0,...,npix()-1].
1526  * @param[in] map Map index [0,...,nmaps()-1].
1527  * @return Flux in pixel.
1528  *
1529  * @exception GException::invalid_value
1530  * No valid sky projection found.
1531  *
1532  * Returns the flux in the pixel with the specified @p index. The flux is
1533  * computed by integrating the intensity over the solid angle subtended by
1534  * the pixel.
1535  * Integration is done by dividing the pixel into wedges, computing the
1536  * intensitites at the cornes of each wedge, averaging these intensities,
1537  * and multiplying it with the solid angle of each wedge. This provides an
1538  * approximation of the true pixel flux which is accurate to better than
1539  * about 5%.
1540  *
1541  * For a HealPix pixelisation, the pixel is divided into 12 wedges. For a
1542  * WCS pixelisation, the pixel is divided into 8 wedges.
1543  *
1544  * @warning
1545  * This method only returns correct values once the skymap is completely
1546  * setup with values. Do not use this method during a setup operation as
1547  * the method uses neighboring pixels for interpolation. If the map is not
1548  * completely setup the neighboring pixels may still be empty, hence the
1549  * flux interpolation will be wrong.
1550  ***************************************************************************/
1551 double GSkyMap::flux(const int& index, const int& map) const
1552 {
1553  // Throw error if WCS is not valid
1554  if (m_proj == NULL) {
1555  std::string msg = "Sky projection has not been defined.";
1556  throw GException::invalid_value(G_FLUX1, msg);
1557  }
1558 
1559  // Determine flux from pixel index.
1560  double flux = this->flux(inx2pix(index), map);
1561 
1562  // Return flux
1563  return flux;
1564 }
1565 
1566 
1567 /***********************************************************************//**
1568  * @brief Returns flux in pixel
1569  *
1570  * @param[in] pixel Sky map pixel.
1571  * @param[in] map Map index [0,...,nmaps()-1].
1572  * @return Flux in pixel (ph/cm2/s).
1573  *
1574  * @exception GException::invalid_value
1575  * No valid sky projection found.
1576  *
1577  * Returns the flux in the specified sky map @p pixel. The flux is computed
1578  * by integrating the intensity over the solid angle subtended by the pixel.
1579  * Integration is done by dividing the pixel into wedges, computing the
1580  * intensitites at the cornes of each wedge, averaging these intensities,
1581  * and multiplying it with the solid angle of each wedge. This provides an
1582  * approximation of the true pixel flux which is accurate to better than
1583  * about 5%.
1584  *
1585  * For a HealPix pixelisation, the pixel is divided into 12 wedges. For a
1586  * WCS pixelisation, the pixel is divided into 8 wedges.
1587  *
1588  * @warning
1589  * This method only returns correct values once the skymap is completely
1590  * setup with values. Do not use this method during a setup operation as
1591  * the method uses neighboring pixels for interpolation. If the map is not
1592  * completely setup the neighboring pixels may still be empty, hence the
1593  * flux interpolation will be wrong.
1594  ***************************************************************************/
1595 double GSkyMap::flux(const GSkyPixel& pixel, const int& map) const
1596 {
1597  // Throw error if WCS is not valid
1598  if (m_proj == NULL) {
1599  std::string msg = "Sky projection has not been defined.";
1600  throw GException::invalid_value(G_FLUX2, msg);
1601  }
1602 
1603  // Initialise flux
1604  double flux = 0.0;
1605 
1606  // Perform flux computation for HealPix map. The pixel is divided into
1607  // 12 wedges and the flux is integrated in each wedge by computing the
1608  // average of the corner intensities multiplied by the solid angle of
1609  // the wedge. This leads to a precision of the order of 1-2% in the
1610  // flux computation.
1611  if (m_proj->size() == 1) {
1612 
1613  // Get pointer on HealPix projection
1614  const GHealpix* healpix = static_cast<const GHealpix*>(projection());
1615 
1616  // Get centre
1617  GSkyDir centre = healpix->pix2dir(pixel);
1618 
1619  // Get boundaries
1620  GSkyDirs boundaries = healpix->boundaries(pixel, 3);
1621 
1622  // Compute intensities
1623  double i0 = this->operator()(centre, map);
1624  double i1 = this->operator()(boundaries[0], map);
1625  double i2 = this->operator()(boundaries[4], map);
1626  double i3 = this->operator()(boundaries[8], map);
1627  double i4 = this->operator()(boundaries[1], map);
1628  double i5 = this->operator()(boundaries[5], map);
1629  double i6 = this->operator()(boundaries[9], map);
1630  double i7 = this->operator()(boundaries[2], map);
1631  double i8 = this->operator()(boundaries[6], map);
1632  double i9 = this->operator()(boundaries[10], map);
1633  double i10 = this->operator()(boundaries[3], map);
1634  double i11 = this->operator()(boundaries[7], map);
1635  double i12 = this->operator()(boundaries[11], map);
1636 
1637  // Compute fluxes in wedges
1638  double flux1 = gammalib::onethird * (i0 + i1 + i2) *
1639  solidangle(centre, boundaries[0], boundaries[4]);
1640  double flux2 = gammalib::onethird * (i0 + i2 + i3) *
1641  solidangle(centre, boundaries[4], boundaries[8]);
1642  double flux3 = gammalib::onethird * (i0 + i3 + i4) *
1643  solidangle(centre, boundaries[8], boundaries[1]);
1644  double flux4 = gammalib::onethird * (i0 + i4 + i5) *
1645  solidangle(centre, boundaries[1], boundaries[5]);
1646  double flux5 = gammalib::onethird * (i0 + i5 + i6) *
1647  solidangle(centre, boundaries[5], boundaries[9]);
1648  double flux6 = gammalib::onethird * (i0 + i6 + i7) *
1649  solidangle(centre, boundaries[9], boundaries[2]);
1650  double flux7 = gammalib::onethird * (i0 + i7 + i8) *
1651  solidangle(centre, boundaries[2], boundaries[6]);
1652  double flux8 = gammalib::onethird * (i0 + i8 + i9) *
1653  solidangle(centre, boundaries[6], boundaries[10]);
1654  double flux9 = gammalib::onethird * (i0 + i9 + i10) *
1655  solidangle(centre, boundaries[10], boundaries[3]);
1656  double flux10 = gammalib::onethird * (i0 + i10 + i11) *
1657  solidangle(centre, boundaries[3], boundaries[7]);
1658  double flux11 = gammalib::onethird * (i0 + i11 + i12) *
1659  solidangle(centre, boundaries[7], boundaries[11]);
1660  double flux12 = gammalib::onethird * (i0 + i12 + i1) *
1661  solidangle(centre, boundaries[11], boundaries[0]);
1662 
1663  // Sum up fluxes
1664  flux = (flux1 + flux2 + flux3 + flux4 +
1665  flux5 + flux6 + flux7 + flux8 +
1666  flux9 + flux10 + flux11 + flux12);
1667 
1668  } // endif: we had a HealPix map
1669 
1670  // ... otherwise perform flux computation for WCS map. The pixel is
1671  // divided into 8 wedges and the flux is integrated in each wedge by
1672  // computing the average of the corner intensities multiplied by the
1673  // solid angle of the wedge. This leads to a precision of the order
1674  // of <1% in the flux computation.
1675  else {
1676 
1677  // Get centre
1678  GSkyDir centre = pix2dir(pixel);
1679 
1680  // Set upper pixel edge just beyond the upper boundary since the
1681  // upper boundary is excluded
1682  const double up = 0.4999999;
1683 
1684  // Get boundaries
1685  GSkyDir boundary1 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()-0.5));
1686  GSkyDir boundary2 = pix2dir(GSkyPixel(pixel.x(), pixel.y()-0.5));
1687  GSkyDir boundary3 = pix2dir(GSkyPixel(pixel.x()+up, pixel.y()-0.5));
1688  GSkyDir boundary4 = pix2dir(GSkyPixel(pixel.x()+up, pixel.y()));
1689  GSkyDir boundary5 = pix2dir(GSkyPixel(pixel.x()+up, pixel.y()+up));
1690  GSkyDir boundary6 = pix2dir(GSkyPixel(pixel.x(), pixel.y()+up));
1691  GSkyDir boundary7 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()+up));
1692  GSkyDir boundary8 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()));
1693 
1694  // Compute intensities
1695  double i0 = this->operator()(centre, map);
1696  double i1 = this->operator()(boundary1, map);
1697  double i2 = this->operator()(boundary2, map);
1698  double i3 = this->operator()(boundary3, map);
1699  double i4 = this->operator()(boundary4, map);
1700  double i5 = this->operator()(boundary5, map);
1701  double i6 = this->operator()(boundary6, map);
1702  double i7 = this->operator()(boundary7, map);
1703  double i8 = this->operator()(boundary8, map);
1704 
1705  // Compute fluxes
1706  double flux1 = gammalib::onethird * (i1 + i2 + i0) *
1707  solidangle(boundary1, boundary2, centre);
1708  double flux2 = gammalib::onethird * (i2 + i3 + i0) *
1709  solidangle(boundary2, boundary3, centre);
1710  double flux3 = gammalib::onethird * (i3 + i4 + i0) *
1711  solidangle(boundary3, boundary4, centre);
1712  double flux4 = gammalib::onethird * (i4 + i5 + i0) *
1713  solidangle(boundary4, boundary5, centre);
1714  double flux5 = gammalib::onethird * (i5 + i6 + i0) *
1715  solidangle(boundary5, boundary6, centre);
1716  double flux6 = gammalib::onethird * (i6 + i7 + i0) *
1717  solidangle(boundary6, boundary7, centre);
1718  double flux7 = gammalib::onethird * (i7 + i8 + i0) *
1719  solidangle(boundary7, boundary8, centre);
1720  double flux8 = gammalib::onethird * (i8 + i1 + i0) *
1721  solidangle(boundary8, boundary1, centre);
1722 
1723  // Sum up fluxes
1724  flux = (flux1 + flux2 + flux3 + flux4 +
1725  flux5 + flux6 + flux7 + flux8);
1726 
1727  } // endelse: we had a WCS map
1728 
1729  // Return flux
1730  return flux;
1731 }
1732 
1733 
1734 /***********************************************************************//**
1735  * @brief Returns flux within sky region
1736  *
1737  * @param[in] region Sky region.
1738  * @param[in] map Map index [0,...,nmaps()-1].
1739  * @return Flux in sky region.
1740  *
1741  * @exception GException::out_of_range
1742  * Map index out of valid range
1743  *
1744  * The method returns the total flux in all sky map pixels that fall within
1745  * the specified region. Containment is tested using the
1746  * GSkyRegion::contains() method. Sky map values are multiplied by the solid
1747  * angle of the pixel.
1748  ***************************************************************************/
1749 double GSkyMap::flux(const GSkyRegion& region, const int& map) const
1750 {
1751  // Throw an error if map index is not in valid range
1752  if (map < 0 || map >= m_num_maps) {
1753  throw GException::out_of_range(G_FLUX3, "Map index",
1754  map, m_num_maps);
1755  }
1756 
1757  // Initialise flux
1758  double flux = 0.0;
1759 
1760  // Loop over all sky pixels and sum flux of all pixels that fall within
1761  // the sky region
1762  for (int i = 0, ipix = m_num_pixels * map; i < m_num_pixels; ++i, ++ipix) {
1763  if (region.contains(inx2dir(i))) {
1764  flux += m_pixels(ipix) * solidangle(i);
1765  }
1766  }
1767 
1768  // Return flux
1769  return flux;
1770 }
1771 
1772 
1773 /***********************************************************************//**
1774  * @brief Returns flux within sky regions
1775  *
1776  * @param[in] regions Sky regions.
1777  * @param[in] map Map index [0,...,nmaps()-1].
1778  * @return Flux in sky region.
1779  *
1780  * @exception GException::out_of_range
1781  * Map index out of valid range
1782  *
1783  * The method returns the total flux in all sky map pixels that fall within
1784  * the specified regions. Containment is tested using the
1785  * GSkyRegions::contains() method. Sky map values are multiplied by the solid
1786  * angle of the pixel.
1787  ***************************************************************************/
1788 double GSkyMap::flux(const GSkyRegions& regions, const int& map) const
1789 {
1790  // Throw an error if map index is not in valid range
1791  if (map < 0 || map >= m_num_maps) {
1792  throw GException::out_of_range(G_FLUX4, "Map index",
1793  map, m_num_maps);
1794  }
1795 
1796  // Initialise flux
1797  double flux = 0.0;
1798 
1799  // Loop over all sky pixels and sum flux of all pixels that fall within
1800  // the sky region
1801  for (int i = 0, ipix = m_num_pixels * map; i < m_num_pixels; ++i, ++ipix) {
1802  if (regions.contains(inx2dir(i))) {
1803  flux += m_pixels(ipix) * solidangle(i);
1804  }
1805  }
1806 
1807  // Return flux
1808  return flux;
1809 }
1810 
1811 
1812 /***********************************************************************//**
1813  * @brief Returns array with total flux for sky maps
1814  *
1815  * @return Array of pixel flux sums.
1816  *
1817  * For a set of n sky maps, the method returns a 1D array with n entries
1818  * that correspond to the sum of the pixel flux in each map.
1819  ***************************************************************************/
1821 {
1822  // Initialise output array
1823  GNdarray flux_array = GNdarray(m_num_maps);
1824 
1825  // Loop over maps and store sums in output array
1826  for (int i = 0, ipix = 0; i < m_num_maps; ++i) {
1827 
1828  // Initialise flux for this map
1829  double flux = 0.0;
1830 
1831  // Compute map flux
1832  for (int k = 0; k < m_num_pixels; ++k, ++ipix) {
1833  double solidangle = this->solidangle(k);
1834  flux += m_pixels(ipix) * solidangle;
1835  }
1836 
1837  // Store map flux
1838  flux_array(i) = flux;
1839 
1840  } // endfor: looped over all maps
1841 
1842  // Return flux array
1843  return flux_array;
1844 }
1845 
1846 
1847 /***********************************************************************//**
1848  * @brief Returns solid angle of pixel
1849  *
1850  * @param[in] index Pixel index [0,...,npix()-1].
1851  * @return Solid angle (steradians)
1852  *
1853  * @exception GException::invalid_value
1854  * No valid sky projection found.
1855  *
1856  * Returns the solid angle of the pixel with the specified @p index.
1857  ***************************************************************************/
1858 double GSkyMap::solidangle(const int& index) const
1859 {
1860  // Throw error if WCS is not valid
1861  if (m_proj == NULL) {
1862  std::string msg = "Sky projection has not been defined.";
1864  }
1865 
1866  // Determine solid angle from pixel index.
1867  double solidangle = this->solidangle(inx2pix(index));
1868 
1869  // Return solid angle
1870  return solidangle;
1871 }
1872 
1873 
1874 /***********************************************************************//**
1875  * @brief Returns solid angle of pixel
1876  *
1877  * @param[in] pixel Sky map pixel.
1878  * @return Solid angle (steradians)
1879  *
1880  * @exception GException::invalid_value
1881  * No valid sky projection found.
1882  *
1883  * Returns the solid angle of the specified sky map @p pixel.
1884  ***************************************************************************/
1885 double GSkyMap::solidangle(const GSkyPixel& pixel) const
1886 {
1887  // Throw error if WCS is not valid
1888  if (m_proj == NULL) {
1889  std::string msg = "Sky projection has not been defined.";
1891  }
1892 
1893  // Initialise solid angle
1894  double solidangle = 0.0;
1895 
1896  // If pixel size matches the projection size then perform a straight
1897  // forward solid angle determination
1898  if (m_proj->size() == pixel.size()) {
1899  solidangle = m_proj->solidangle(pixel);
1900  }
1901 
1902  // ... otherwise, if we have a 2D projection but a 1D pixel then
1903  // interpret the pixel as the linear index in the pixel array
1904  else if (m_proj->size() == 2) {
1905  solidangle = m_proj->solidangle(GSkyPixel(inx2pix(int(pixel))));
1906  }
1907 
1908  // ... otherwise we have a 1D projection but a 2D pixel. There is
1909  // no unambiguous way to handle this case, hence we throw an exception
1910  else {
1911  std::string msg = "A 2-dimensional sky map pixel "+pixel.print()+
1912  " is used to determine the solid angle for"
1913  " the 1-dimensional sky projection \""+
1914  m_proj->name()+"\"\n"
1915  "Please specify a 1-dimensional sky map pixel.";
1917  }
1918 
1919  // Return solid angle
1920  return solidangle;
1921 }
1922 
1923 
1924 /***********************************************************************//**
1925  * @brief Returns solid angle within sky region
1926  *
1927  * @param[in] region Sky region.
1928  * @return Solid angle within sky region.
1929  *
1930  * The method returns the total solid angle of all sky map pixels that fall
1931  * within the specified region. Containment is tested using the
1932  * GSkyRegion::contains() method.
1933  ***************************************************************************/
1934 double GSkyMap::solidangle(const GSkyRegion& region) const
1935 {
1936  // Initialise flux
1937  double solidangle = 0.0;
1938 
1939  // Loop over all sky pixels and sum solid angle of all pixels that fall
1940  // within the sky region
1941  for (int i = 0; i < m_num_pixels; ++i) {
1942  if (region.contains(inx2dir(i))) {
1943  solidangle += this->solidangle(i);
1944  }
1945  }
1946 
1947  // Return solidangle
1948  return solidangle;
1949 }
1950 
1951 
1952 /***********************************************************************//**
1953  * @brief Returns solid angle within sky regions
1954  *
1955  * @param[in] regions Sky regions.
1956  * @return Solid angle within sky region.
1957  *
1958  * The method returns the total solid angle of all sky map pixels that fall
1959  * within the specified regions. Containment is tested using the
1960  * GSkyRegions::contains() method.
1961  ***************************************************************************/
1962 double GSkyMap::solidangle(const GSkyRegions& regions) const
1963 {
1964  // Initialise flux
1965  double solidangle = 0.0;
1966 
1967  // Loop over all sky pixels and sum solid angle of all pixels that fall
1968  // within the sky regions
1969  for (int i = 0; i < m_num_pixels; ++i) {
1970  if (regions.contains(inx2dir(i))) {
1971  solidangle += this->solidangle(i);
1972  }
1973  }
1974 
1975  // Return solidangle
1976  return solidangle;
1977 }
1978 
1979 
1980 /***********************************************************************//**
1981  * @brief Set sky projection
1982  *
1983  * @param[in] proj Sky projection.
1984  *
1985  * Sets the projection from celestial to pixel coordinates. The method
1986  * performs a deep copy of @p proj, allowing to destroy the argument after
1987  * using the method.
1988  *
1989  * Warning: this method may corrupt the GSkyMap object as it allows assigning
1990  * for example a 1D projection to a 2D skymap. Please use this method only
1991  * when you know what you're doing.
1992  *
1993  * @todo We may restrict this method to not allow changing the projection
1994  * dimension.
1995  ***************************************************************************/
1997 {
1998  // Free existing projection only if it differs from current projection.
1999  // This prevents unintential deallocation of the argument
2000  if ((m_proj != NULL) && (m_proj != &proj)) {
2001  delete m_proj;
2002  }
2003 
2004  // Clone input WCS
2005  m_proj = proj.clone();
2006 
2007  // Return
2008  return;
2009 }
2010 
2011 
2012 /***********************************************************************//**
2013  * @brief Checks if sky direction falls in map
2014  *
2015  * @param[in] dir Sky direction.
2016  * @return True if sky direction falls into map, false otherwise.
2017  *
2018  * This method checks if the specified sky direction falls within the pixels
2019  * covered by the skymap. The method uses the dir2xy method to convert the
2020  * sky direction into 2D pixel indices, and then checks whether the pixel
2021  * indices fall in the skymap.
2022  ***************************************************************************/
2023 bool GSkyMap::contains(const GSkyDir& dir) const
2024 {
2025  // Initialise containment
2026  bool contains = false;
2027 
2028  // Convert sky direction into sky pixel
2029  try {
2030 
2031  // Convert sky direction into sky pixel
2032  GSkyPixel pixel = dir2pix(dir);
2033 
2034  // Check pixel containment
2035  contains = this->contains(pixel);
2036 
2037  }
2039  contains = false;
2040  }
2041 
2042  // Return containment flag
2043  return contains;
2044 }
2045 
2046 
2047 /***********************************************************************//**
2048  * @brief Checks if sky map pixel falls in map
2049  *
2050  * @param[in] pixel Sky map pixel.
2051  * @return Trus if pixels is within map, false otherwise.
2052  *
2053  * Checks whether the specified sky map @p pixel falls within the skymap
2054  * or not. A pixel is considered to fall in the sky may if the value is
2055  * contained in the interval
2056  *
2057  * \f[
2058  * [-0.5, n-0.5[
2059  * \f]
2060  *
2061  * where \f$n\f$ is the number of sky map pixels. Note that the upper bound
2062  * is excluded from the interval.
2063  ***************************************************************************/
2064 bool GSkyMap::contains(const GSkyPixel& pixel) const
2065 {
2066  // Initialise containment flag
2067  bool inmap = false;
2068 
2069  // Test 1D pixels
2070  if (pixel.is_1D()) {
2071  if (pixel.index() >= -0.5 &&
2072  pixel.index() < double(m_num_pixels)-0.5) {
2073  inmap = true;
2074  }
2075  }
2076 
2077  // Test 2D pixels
2078  else if (pixel.is_2D()) {
2079 
2080  // If pixel is in range then set containment flag to true
2081  if ((pixel.x() >= -0.5 && pixel.x() < double(m_num_x)-0.5) &&
2082  (pixel.y() >= -0.5 && pixel.y() < double(m_num_y)-0.5)) {
2083  inmap = true;
2084  }
2085 
2086  }
2087 
2088  // Return containment flag
2089  return inmap;
2090 }
2091 
2092 
2093 /***********************************************************************//**
2094  * @brief Checks whether a region overlaps with this map
2095  *
2096  * @param[in] region Region
2097  * @return True if region overlaps with map, false otherwise
2098  *
2099  * @exception GException::feature_not_implemented
2100  * Region is not a circular sky region
2101  ***************************************************************************/
2102 bool GSkyMap::overlaps(const GSkyRegion& region) const
2103 {
2104  // Initialise overlap
2105  bool overlap = false;
2106 
2107  // If the region is a circular region than test on overlap
2108  if (region.type() == "Circle") {
2109  const GSkyRegionCircle* circle =
2110  static_cast<const GSkyRegionCircle*>(&region);
2111  overlap = overlaps_circle(*circle);
2112  }
2113 
2114  // ... otherwise signal that method is not implemented
2115  else {
2116  std::string msg = "Method not implemented for non-circular regions.";
2118  }
2119 
2120  // Return overlap flag
2121  return overlap;
2122 }
2123 
2124 
2125 /***********************************************************************//**
2126  * @brief Extract maps into a new sky map object
2127  *
2128  * @param[in] map First map to extract
2129  * @param[in] nmaps Number of maps to extract
2130  * @return Extracted map(s).
2131  *
2132  * @exception GException::out_of_range
2133  * First map index outside valid range
2134  * @exception GException::invalid_argument
2135  * Requested number of maps are not available.
2136  *
2137  * Extracts @p nmaps sky maps starting from @p map from the sky map.
2138  ***************************************************************************/
2139 GSkyMap GSkyMap::extract(const int& map, const int& nmaps) const
2140 {
2141  // Throw an exception if the first map index is invalid
2142  if (map < 0 || map >= m_num_maps) {
2143  throw GException::out_of_range(G_EXTRACT, "Sky map index", map,
2144  m_num_maps);
2145  }
2146 
2147  // Throw an exception if the number of maps is invalid
2148  if (nmaps < 0) {
2149  std::string msg = "The number of maps to extract cannot be negative "
2150  "(nmaps="+gammalib::str(nmaps)+").";
2152  }
2153  if (nmaps > m_num_maps-map) {
2154  std::string msg = "The number of maps to extract ("+
2155  gammalib::str(nmaps)+") exceeds the number of maps "
2156  "that are available ("+
2157  gammalib::str(m_num_maps-map)+").";
2159  }
2160 
2161  // Allocate extracted pixels
2162  GNdarray extracted_pixels(m_num_pixels, nmaps);
2163 
2164  // Extract pixels
2165  int n_size = m_num_pixels * nmaps;
2166  const double *src = m_pixels.data() + map*m_num_pixels;
2167  double *dst = extracted_pixels.data();
2168  for (int i = 0; i < n_size; ++i) {
2169  *dst++ = *src++;
2170  }
2171 
2172  // Create a copy of the map
2173  GSkyMap result = *this;
2174 
2175  // Attach copied pixels to the map
2176  result.m_pixels = extracted_pixels;
2177 
2178  // Set number of maps
2179  result.m_num_maps = nmaps;
2180 
2181  // Reset shape of maps since we do not know how to arrange the maps
2182  result.m_shape.clear();
2183  result.m_shape.push_back(nmaps);
2184 
2185  // Return map
2186  return result;
2187 }
2188 
2189 
2190 /***********************************************************************//**
2191  * @brief Extract a sub-range of pixels in the map in x,y
2192  *
2193  * @param[in] startx Starting bin in X (inclusive)
2194  * @param[in] stopx Last bin in X (inclusive)
2195  * @param[in] starty Starting bin in Y (inclusive)
2196  * @param[in] stopy Last bin in Y (inclusive)
2197  * @return Skymap that overlaps with supplied exclusions
2198  *
2199  * @exception GException::invalid_argument
2200  * Method not valid for HPX projection
2201  *
2202  * This method creates a new skymap consisting of all pixels in the map in the
2203  * range [startx,stopx] and [starty,stopy]. The boundary values provided are
2204  * inclusive and the actual projection is unchanged. The number of maps in the
2205  * skymap is also unchanged.
2206  *
2207  * NOTE: The pixel ranges are indexed from 0, as in
2208  ***************************************************************************/
2209 GSkyMap GSkyMap::extract(const int& startx, const int& stopx,
2210  const int& starty, const int& stopy) const
2211 {
2212  // Make sure this isn't a 'HealPix' map
2213  if (m_proj->code() == "HPX") {
2214  std::string msg = "Method not valid for \"HPX\" projection. Please "
2215  "specify WCS projection.";
2217  }
2218  else if ((startx > stopx) || (starty > stopy)) {
2220  "Invalid x,y range provided");
2221  }
2222 
2223  // Define the actual range of pixels
2224  int first_x = (startx < 0) ? 0 : startx;
2225  int last_x = (stopx >= m_num_x) ? m_num_x-1 : stopx;
2226  int first_y = (starty < 0) ? 0 : starty;
2227  int last_y = (stopy >= m_num_y) ? m_num_y-1 : stopy;
2228  int nxpix = last_x - first_x + 1;
2229  int nypix = last_y - first_y + 1;
2230  int npixels = nxpix * nypix;
2231 
2232  // Create an array of values to be plotted
2233  GNdarray new_pixels(nxpix, nypix, m_num_maps);
2234 
2235  // Get new pixel values
2236  for (int ix = 0; ix < nxpix; ++ix) {
2237  int ix_src = ix + first_x;
2238  for (int iy = 0; iy < nypix; ++iy) {
2239  int iy_src = iy + first_y;
2240  for (int imap = 0; imap < m_num_maps; ++imap) {
2241  new_pixels(ix, iy, imap) = m_pixels(ix_src, iy_src, imap);
2242  }
2243  }
2244  }
2245 
2246  // Create a copy of this map
2247  GSkyMap result = *this;
2248 
2249  // Update the projection
2250  GWcs* proj = static_cast<GWcs*>(result.m_proj);
2251  proj->set(proj->coordsys(),
2252  proj->crval(0), proj->crval(1),
2253  proj->crpix(0) - first_x, // Offset reference pixel
2254  proj->crpix(1) - first_y, // Offset reference pixel
2255  proj->cdelt(0), proj->cdelt(1));
2256 
2257  // Update the pixel information
2258  result.m_pixels = new_pixels;
2259  result.m_num_x = nxpix;
2260  result.m_num_y = nypix;
2261  result.m_num_pixels = npixels;
2262 
2263  // Re-initialise computation cache
2264  result.m_hascache = false;
2265  result.m_contained = false;
2266  result.m_last_dir.clear();
2267  result.m_interpol.clear();
2268 
2269  // Return map
2270  return result;
2271 }
2272 
2273 
2274 /***********************************************************************//**
2275  * @brief Extract the spatial portion of the maps that overlap @p inclusions
2276  *
2277  * @param[in] inclusions List of GSkyRegion objects
2278  * @return Skymap that overlaps with supplied exclusions
2279  *
2280  * @exception GException::invalid_argument
2281  * Method not valid for HPX projection
2282  *
2283  * This method computes the x,y range that covers the regions in @p inclusions,
2284  * then returns that range as a new GSkyMap.
2285  *
2286  * NOTE: A 1-pixel border is added to the returned map to provide better
2287  * coverage in the returned map.
2288  ***************************************************************************/
2289 GSkyMap GSkyMap::extract(const GSkyRegions& inclusions) const
2290 {
2291  // Make sure this isn't a 'HealPix' map
2292  if (m_proj->code() == "HPX") {
2293  std::string msg = "Method not valid for \"HPX\" projection. Please "
2294  "specify WCS projection.";
2296  }
2297 
2298  // Preliminary variables
2299  int startx = nx()+1;
2300  int stopx = 0;
2301  int starty = ny()+1;
2302  int stopy = 0;
2303 
2304  // Loop on pixels
2305  for (int i = 0; i < npix(); ++i) {
2306 
2307  // Get the pixel and direction
2308  GSkyPixel pixel = inx2pix(i);
2309  GSkyDir pixdir = pix2dir(pixel);
2310 
2311  // Check if the pixel is covered by 'inclusions'
2312  if (inclusions.contains(pixdir)) {
2313 
2314  // Get the x,y pixel
2315  int pix_x = int(pixel.x());
2316  int pix_y = int(pixel.y());
2317 
2318  // Update x,y pixel ranges
2319  startx = (pix_x < startx) ? pix_x : startx;
2320  stopx = (pix_x > stopx) ? pix_x : stopx;
2321  starty = (pix_y < starty) ? pix_y : starty;
2322  stopy = (pix_y > stopy) ? pix_y : stopy;
2323 
2324  } // endif: pixel was covered by inclusions
2325 
2326  } // endfor: looped over all pixels
2327 
2328  // Trim the skymap with 1 pixel buffer all around.
2329  return (this->extract(startx-1, stopx+1, starty-1, stopy+1));
2330 }
2331 
2332 
2333 /***********************************************************************//**
2334  * @brief Stack all maps into a single map
2335  *
2336  * The methods replaces the sky map by a version with only a single map
2337  * by summing over the pixel values for all maps. If the sky map has no
2338  * pixels or there is only a single map in the object, the method does
2339  * nothing.
2340  ***************************************************************************/
2342 {
2343  // Continue only if the map has pixels and if there is more than one map
2344  if (m_num_pixels > 0 && m_num_maps > 1) {
2345 
2346  // Initialise stacked pixels
2347  GNdarray stacked_pixels;
2348 
2349  // Allocate stacked pixels
2350  if (m_proj->code() == "HPX") {
2351  stacked_pixels = GNdarray(m_num_pixels);
2352  }
2353  else {
2354  stacked_pixels = GNdarray(m_num_x, m_num_y);
2355  }
2356 
2357  // Stack map and save in pixels
2358  for (int i = 0; i < m_num_pixels; ++i) {
2359  double sum = 0.0;
2360  for (int imap = 0; imap < m_num_maps; ++imap) {
2361  sum += (*this)(i,imap);
2362  }
2363  stacked_pixels(i) = sum;
2364  }
2365 
2366  // Set pixels to stacked pixels
2367  m_pixels = stacked_pixels;
2368 
2369  // Set number of maps to 1
2370  m_num_maps = 1;
2371 
2372  // Reset shape of maps
2373  m_shape.clear();
2374  m_shape.push_back(m_num_maps);
2375 
2376  } // endif: map had pixels
2377 
2378  // Return
2379  return;
2380 }
2381 
2382 
2383 /***********************************************************************//**
2384  * @brief Return sky region circle that encloses the sky map
2385  *
2386  * @return Enclosing sky region circle.
2387  *
2388  * Returns a sky region circle that encloses the sky map.
2389  *
2390  * For a sky map in HealPix projection the method returns a region circle
2391  * that encloses the full sky, with a centre at Right Ascension and
2392  * Declination of 0 and a radius of 180 degrees.
2393  *
2394  * For a sky map in World Coordinate projection the method returns a region
2395  * circle that is centred on the central pixel of the sky map (i.e. the
2396  * pixel with indices \f$(n_x/2,n_y/2)\f$, where \f$n_x\f$ and \f$n_y\f$ are
2397  * the number of pixels in the X- and Y-drection) and a radius that is
2398  * determined from the maximum radial distance of all pixels with respect
2399  * to the central pixel.
2400  ***************************************************************************/
2402 {
2403  // Initialise sky region circle
2404  GSkyRegionCircle region;
2405 
2406  // If the sky map is in HealPix projection then set radius to 180 deg
2407  if (projection()->code() == "HPX") {
2408  region = GSkyRegionCircle(0.0, 0.0, 180.0);
2409  }
2410 
2411  // ... otherwise sky map is in World Coordinate Projection
2412  else {
2413 
2414  // Initialise maximum radius
2415  double max_radius = 0.0;
2416  int max_index = -1;
2417 
2418  // Determine map centre
2419  GSkyPixel pixel(nx()/2.0, ny()/2.0);
2420  GSkyDir centre = pix2dir(pixel);
2421 
2422  // Determine maximum pixel distance with respect to map centre
2423  for (int i = 0; i < npix(); ++i) {
2424  double radius = inx2dir(i).dist_deg(centre);
2425  if (radius > max_radius) {
2426  max_radius = radius;
2427  max_index = i;
2428  }
2429  }
2430 
2431  // Consider edges of pixel with maximum distance
2432  if (max_index != -1) {
2433 
2434  // Set upper pixel edge just beyond the upper boundary since the
2435  // upper boundary is excluded
2436  const double up = 0.4999999;
2437 
2438  // Get pixel with maximum distance
2439  GSkyPixel pixel = inx2pix(max_index);
2440 
2441  // Get boundaries of pixel
2442  GSkyDir boundary1 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()-0.5));
2443  GSkyDir boundary2 = pix2dir(GSkyPixel(pixel.x(), pixel.y()-0.5));
2444  GSkyDir boundary3 = pix2dir(GSkyPixel(pixel.x()+up, pixel.y()-0.5));
2445  GSkyDir boundary4 = pix2dir(GSkyPixel(pixel.x()+up, pixel.y()));
2446  GSkyDir boundary5 = pix2dir(GSkyPixel(pixel.x()+up, pixel.y()+up));
2447  GSkyDir boundary6 = pix2dir(GSkyPixel(pixel.x(), pixel.y()+up));
2448  GSkyDir boundary7 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()+up));
2449  GSkyDir boundary8 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()));
2450 
2451  // Get radii
2452  double radius1 = boundary1.dist_deg(centre);
2453  double radius2 = boundary2.dist_deg(centre);
2454  double radius3 = boundary3.dist_deg(centre);
2455  double radius4 = boundary4.dist_deg(centre);
2456  double radius5 = boundary5.dist_deg(centre);
2457  double radius6 = boundary6.dist_deg(centre);
2458  double radius7 = boundary7.dist_deg(centre);
2459  double radius8 = boundary8.dist_deg(centre);
2460 
2461  // Determine maximum radius
2462  if (radius1 > max_radius) {
2463  max_radius = radius1;
2464  }
2465  if (radius2 > max_radius) {
2466  max_radius = radius2;
2467  }
2468  if (radius3 > max_radius) {
2469  max_radius = radius3;
2470  }
2471  if (radius4 > max_radius) {
2472  max_radius = radius4;
2473  }
2474  if (radius5 > max_radius) {
2475  max_radius = radius5;
2476  }
2477  if (radius6 > max_radius) {
2478  max_radius = radius6;
2479  }
2480  if (radius7 > max_radius) {
2481  max_radius = radius7;
2482  }
2483  if (radius8 > max_radius) {
2484  max_radius = radius8;
2485  }
2486 
2487  } // endif: considered edges of pixel
2488 
2489  // Set sky region
2490  region = GSkyRegionCircle(centre, max_radius);
2491 
2492  } // endelse: sky map was in World Coordinate Projection
2493 
2494  // Return region
2495  return region;
2496 }
2497 
2498 
2499 /***********************************************************************//**
2500  * @brief Load skymap from FITS file.
2501  *
2502  * @param[in] filename FITS file name..
2503  *
2504  * Loads HEALPix and non HEALPix skymaps. First searches for HEALPix map in
2505  * FITS file by scanning all HDUs for PIXTYPE=HEALPIX. If no HEALPix map has
2506  * been found then search load first non-empty image.
2507  *
2508  * @todo Do we have to restrict a HEALPix map to a BinTable and a WCS map
2509  * to a Double precision image???
2510  ***************************************************************************/
2511 void GSkyMap::load(const GFilename& filename)
2512 {
2513  // Clear sky map
2514  clear();
2515 
2516  // Initialize load flag
2517  bool loaded(false);
2518 
2519  // Open FITS file
2520  GFits fits(filename);
2521 
2522  // If an extension name is specified then first try loading that
2523  // extension
2524  if (filename.has_extname()) {
2525  const GFitsHDU& hdu = *fits.at(filename.extname());
2526  if (is_healpix(hdu)) {
2527  read_healpix(static_cast<const GFitsTable&>(hdu));
2528  loaded = true;
2529  }
2530  else if (is_wcs(hdu)) {
2531  read_wcs(static_cast<const GFitsImage&>(hdu));
2532  loaded = true;
2533  }
2534  }
2535 
2536  // ... otherwise is an extension number is specified then try loading
2537  // the corresponding HDU
2538  else if (filename.has_extno()) {
2539  const GFitsHDU& hdu = *fits.at(filename.extno());
2540  if (is_healpix(hdu)) {
2541  read_healpix(static_cast<const GFitsTable&>(hdu));
2542  loaded = true;
2543  }
2544  else if (is_wcs(hdu)) {
2545  read_wcs(static_cast<const GFitsImage&>(hdu));
2546  loaded = true;
2547  }
2548  }
2549 
2550  // If no map has yet been loaded then scan the file for an appropriate
2551  // HDU and read sky map from the first suitable HDU
2552  if (!loaded) {
2553 
2554  // Get number of HDUs
2555  int num = fits.size();
2556 
2557  // First search for HEALPix extension. We can skip the first
2558  // extension since this is always an image and a HEALPix map
2559  // is stored in a binary table
2560  for (int extno = 1; extno < num; ++extno) {
2561 
2562  // Get reference to HDU
2563  const GFitsHDU& hdu = *fits.at(extno);
2564 
2565  // If HDU is HEALPix then read data
2566  if (is_healpix(hdu)) {
2567  read_healpix(static_cast<const GFitsTable&>(hdu));
2568  loaded = true;
2569  break;
2570  }
2571 
2572  } // endfor: looped over HDUs
2573 
2574  // If we have not found a HEALPIX map then search now for an
2575  // image.
2576  if (!loaded) {
2577  for (int extno = 0; extno < num; ++extno) {
2578 
2579  // Get referene to HDU
2580  const GFitsHDU& hdu = *fits.at(extno);
2581 
2582  // If HDU is WCS then read data
2583  if (is_wcs(hdu)) {
2584  read_wcs(static_cast<const GFitsImage&>(hdu));
2585  //loaded = true;
2586  break;
2587  }
2588 
2589  } // endfor: looped over HDUs
2590  } // endif: no sky map yet loaded
2591  } // endif: no sky map yet loaded
2592 
2593  // Close FITS file
2594  fits.close();
2595 
2596  // Return
2597  return;
2598 }
2599 
2600 
2601 /***********************************************************************//**
2602  * @brief Save sky map into FITS file
2603  *
2604  * @param[in] filename FITS file name.
2605  * @param[in] clobber Overwrite existing file?
2606  *
2607  * Saves the sky map into a FITS file. If the file exists already and the
2608  * @p clobber parameter is true, the method will overwrite the content of
2609  * the existing file. Otherwise, an exception is thrown.
2610  ***************************************************************************/
2611 void GSkyMap::save(const GFilename& filename, const bool& clobber) const
2612 {
2613  // Continue only if we have data to save
2614  if (m_proj != NULL) {
2615 
2616  // Initialise HDU pointer
2617  GFitsHDU* hdu = NULL;
2618 
2619  // Case A: Skymap is Healpix
2620  if (m_proj->code() == "HPX") {
2621  hdu = create_healpix_hdu();
2622  }
2623 
2624  // Case B: Skymap is not Healpix
2625  else {
2626  hdu = create_wcs_hdu();
2627  }
2628 
2629  // If we have a valid HDU then save it now to the FITS file
2630  if (hdu != NULL) {
2631 
2632  // Set extension name
2633  if (filename.has_extname()) {
2634  hdu->extname(filename.extname());
2635  }
2636 
2637  // Initialise empty FITS file
2638  GFits fits;
2639 
2640  // Append sky map to FITS file
2641  fits.append(*hdu);
2642 
2643  // Save FITS file (without extension name which was extracted
2644  // earlier and set in the HDU)
2645  fits.saveto(filename.url(), clobber);
2646 
2647  // Delete HDU
2648  delete hdu;
2649 
2650  } // endif: HDU was valid
2651 
2652  } // endif: we had data to save
2653 
2654  // Return
2655  return;
2656 }
2657 
2658 
2659 /***********************************************************************//**
2660  * @brief Read skymap from FITS HDU
2661  *
2662  * @param[in] hdu FITS HDU.
2663  ***************************************************************************/
2664 void GSkyMap::read(const GFitsHDU& hdu)
2665 {
2666  // Free memory and initialise members
2667  free_members();
2668  init_members();
2669 
2670  // If HDU is a Healpix HDU then load a Healpix map
2671  if (is_healpix(hdu)) {
2672  read_healpix(static_cast<const GFitsTable&>(hdu));
2673  }
2674 
2675  // ... otherwise try loading as WCS map
2676  else if (is_wcs(hdu)) {
2677  read_wcs(static_cast<const GFitsImage&>(hdu));
2678  }
2679 
2680  // Return
2681  return;
2682 }
2683 
2684 
2685 /***********************************************************************//**
2686  * @brief Write sky map into FITS file
2687  *
2688  * @param[in] file FITS file pointer.
2689  * @param[in] extname Sky map extension name.
2690  * @return Pointer to written HDU.
2691  *
2692  * Write the sky map into a FITS file. Optionally, the extension name of
2693  * the FITS HDU can be specified using the @p extname parameter. The method
2694  * returns a pointer to the appended HDU. If no HDU has been appended, for
2695  * example because the sky map is empty, the method will return NULL.
2696  ***************************************************************************/
2697 GFitsHDU* GSkyMap::write(GFits& file, const std::string& extname) const
2698 {
2699  // Initialise pointer to appended HDU
2700  GFitsHDU* hdu_appended = NULL;
2701 
2702  // Continue only if we have data to save
2703  if (m_proj != NULL) {
2704 
2705  // Initialise local HDU pointer
2706  GFitsHDU* hdu = NULL;
2707 
2708  // Case A: Skymap is Healpix
2709  if (m_proj->code() == "HPX") {
2710  hdu = create_healpix_hdu();
2711  }
2712 
2713  // Case B: Skymap is not Healpix
2714  else {
2715  hdu = create_wcs_hdu();
2716  }
2717 
2718  // Append HDU to FITS file.
2719  if (hdu != NULL) {
2720 
2721  // Optionally set extension name
2722  if (extname.length() > 0) {
2723  hdu->extname(extname);
2724  }
2725 
2726  // Append HDU and recover the pointer to the appended HDU
2727  hdu_appended = file.append(*hdu);
2728 
2729  // Delete HDU
2730  delete hdu;
2731 
2732  } // endif: there was a valid HDU
2733 
2734  } // endif: we had data to save
2735 
2736  // Return pointer to appended HDU
2737  return hdu_appended;
2738 }
2739 
2740 
2741 /***********************************************************************//**
2742  * @brief Publish sky map
2743  *
2744  * @param[in] name Name of sky map.
2745  *
2746  * Publishes the sky map on a Virtual Observatory Hub. If no Hub is currently
2747  * active, the method will start a new Hub. If on sky map has been allocated
2748  * the method does nothing.
2749  ***************************************************************************/
2750 void GSkyMap::publish(const std::string& name) const
2751 {
2752  // Create FITS file containing the sky map
2753  GFits fits;
2754 
2755  // Write sky map into FITS file and return a pointer to the written HDU.
2756  // This method returns a NULL pointer if no sky map projection exists,
2757  // which typically is true for empty sky maps.
2758  GFitsHDU* hdu = write(fits);
2759 
2760  // Continue only if HDU pointer is valid
2761  if (hdu != NULL) {
2762 
2763  // Optionally set extension name
2764  if (!name.empty()) {
2765  hdu->extname(name);
2766  }
2767 
2768  // Create VO Client
2769  GVOClient client;
2770 
2771  // Publish map using VO client
2772  client.publish(*hdu);
2773 
2774  } // endif: HDU pointer was valid
2775 
2776  // Return
2777  return;
2778 }
2779 
2780 
2781 /***********************************************************************//**
2782  * @brief Print sky map
2783  *
2784  * @param[in] chatter Chattiness.
2785  * @return String containing sky map information.
2786  ***************************************************************************/
2787 std::string GSkyMap::print(const GChatter& chatter) const
2788 {
2789  // Initialise result string
2790  std::string result;
2791 
2792  // Continue only if chatter is not silent
2793  if (chatter != SILENT) {
2794 
2795  // Append header
2796  result.append("=== GSkyMap ===");
2797 
2798  // Append number of pixels and number of maps
2799  result.append("\n"+gammalib::parformat("Number of pixels"));
2800  result.append(gammalib::str(m_num_pixels));
2801 
2802  // Append WCS dimension information
2803  if (m_proj != NULL && m_proj->size() == 2) {
2804  result.append("\n"+gammalib::parformat("X axis dimension"));
2805  result.append(gammalib::str(m_num_x));
2806  result.append("\n"+gammalib::parformat("Y axis dimension"));
2807  result.append(gammalib::str(m_num_y));
2808  }
2809 
2810  // Append number of maps
2811  result.append("\n"+gammalib::parformat("Number of maps"));
2812  result.append(gammalib::str(m_num_maps));
2813 
2814  // Append shape of maps
2815  result.append("\n"+gammalib::parformat("Shape of maps"));
2816  if (m_shape.size() > 0) {
2817  std::string shape = "(";
2818  for (int i = 0; i < m_shape.size(); ++i) {
2819  if (i > 0) {
2820  shape += ", ";
2821  }
2822  shape += gammalib::str(m_shape[i]);
2823  }
2824  shape += ")";
2825  result.append(shape);
2826  }
2827  else {
2828  result.append("not defined");
2829  }
2830 
2831  // Append sky projection information
2832  if (gammalib::reduce(chatter) > SILENT) {
2833  if (m_proj != NULL) {
2834  result.append("\n"+m_proj->print(gammalib::reduce(chatter)));
2835  }
2836  else {
2837  result.append("\n"+gammalib::parformat("Sky projection"));
2838  result.append("not defined");
2839  }
2840  }
2841 
2842  } // endif: chatter was not silent
2843 
2844  // Return result
2845  return result;
2846 }
2847 
2848 /*==========================================================================
2849  = =
2850  = Private methods =
2851  = =
2852  ==========================================================================*/
2853 
2854 /***********************************************************************//**
2855  * @brief Initialise class members
2856  ***************************************************************************/
2858 {
2859  // Initialise members
2860  m_num_pixels = 0;
2861  m_num_maps = 0;
2862  m_num_x = 0;
2863  m_num_y = 0;
2864  m_shape.clear();
2865  m_proj = NULL;
2866  m_pixels.clear();
2867 
2868  // Initialise computation cache
2869  m_hascache = false;
2870  m_contained = false;
2871  m_last_dir.clear();
2872  m_interpol.clear();
2873 
2874  // Return
2875  return;
2876 }
2877 
2878 
2879 /***********************************************************************//**
2880  * @brief Copy class members
2881  *
2882  * @param[in] map Sky map.
2883  ***************************************************************************/
2885 {
2886  // Copy attributes
2887  m_num_pixels = map.m_num_pixels;
2888  m_num_maps = map.m_num_maps;
2889  m_num_x = map.m_num_x;
2890  m_num_y = map.m_num_y;
2891  m_shape = map.m_shape;
2892  m_pixels = map.m_pixels;
2893 
2894  // Copy computation cache
2895  m_hascache = map.m_hascache;
2896  m_contained = map.m_contained;
2897  m_last_dir = map.m_last_dir;
2898  m_interpol = map.m_interpol;
2899 
2900  // Clone sky projection if it is valid
2901  if (map.m_proj != NULL) m_proj = map.m_proj->clone();
2902 
2903  // Return
2904  return;
2905 }
2906 
2907 
2908 /***********************************************************************//**
2909  * @brief Delete class members
2910  ***************************************************************************/
2912 {
2913  // Free memory
2914  if (m_proj != NULL) {
2915  delete m_proj;
2916  }
2917 
2918  // Signal free pointers
2919  m_proj = NULL;
2920 
2921  // Return
2922  return;
2923 }
2924 
2925 
2926 /***********************************************************************//**
2927  * @brief Set World Coordinate System
2928  *
2929  * @param[in] wcs World Coordinate System code.
2930  * @param[in] coords Coordinate system.
2931  * @param[in] crval1 X value of reference pixel.
2932  * @param[in] crval2 Y value of reference pixel.
2933  * @param[in] crpix1 X index of reference pixel.
2934  * @param[in] crpix2 Y index of reference pixel.
2935  * @param[in] cdelt1 Increment in x direction at reference pixel (deg).
2936  * @param[in] cdelt2 Increment in y direction at reference pixel (deg).
2937  *
2938  * @exception GException::invalid_argument
2939  * Invalid wcs parameter (World Coordinate System not known).
2940  *
2941  * This method sets the WCS projection pointer based on the WCS code and
2942  * sky map parameters. It makes use of the GWcsRegistry class to allocate
2943  * the correct derived class. Note that this method does not support the
2944  * HPX projection.
2945  ***************************************************************************/
2946 void GSkyMap::set_wcs(const std::string& wcs,
2947  const std::string& coords,
2948  const double& crval1,
2949  const double& crval2,
2950  const double& crpix1,
2951  const double& crpix2,
2952  const double& cdelt1,
2953  const double& cdelt2)
2954 {
2955  // Convert WCS to upper case
2956  std::string uwcs = gammalib::toupper(wcs);
2957 
2958  // Check if HPX was requested (since this is not allowed)
2959  if (uwcs == "HPX") {
2960  std::string msg = "Method not valid for \"HPX\" projection. Please "
2961  "specify WCS projections.";
2963  }
2964 
2965  // ... otherwise get projection from registry
2966  else {
2967  // Allocate WCS registry
2968  GWcsRegistry registry;
2969 
2970  // Allocate projection from registry
2971  GWcs* projection = registry.alloc(uwcs);
2972  m_proj = projection;
2973 
2974  // Signal if projection type is not known
2975  if (projection == NULL) {
2976  std::string msg = "WCS projection code \""+uwcs+"\" not known. "
2977  "Please specify one of the following WCS "
2978  "projections: "+registry.content()+".";
2980  }
2981 
2982  // Setup WCS
2983  projection->set(coords, crval1, crval2, crpix1, crpix2,
2984  cdelt1, cdelt2);
2985 
2986  } // endelse: got projection from registry
2987 
2988  // Return
2989  return;
2990 }
2991 
2992 
2993 /***********************************************************************//**
2994  * @brief Read Healpix data from FITS table.
2995  *
2996  * @param[in] table FITS table.
2997  *
2998  * @exception GException::runtime_error
2999  * Incompatible HealPix projection encountered.
3000  *
3001  * HEALPix data may be stored in various formats depending on the
3002  * application that has writted the data. HEALPix IDL, for example, may
3003  * store the data in vectors of length 1024 if the number of pixels is
3004  * a multiple of 1024. On the other hand, vectors may also be used to store
3005  * several HEALPix maps into a single column. Alternatively, multiple maps
3006  * may be stored in multiple columns.
3007  ***************************************************************************/
3009 {
3010  // Determine number of rows and columns in table
3011  int nrows = table.nrows();
3012  int ncols = table.ncols();
3013  #if defined(G_READ_HEALPIX_DEBUG)
3014  std::cout << "nrows=" << nrows << " ncols=" << ncols << std::endl;
3015  #endif
3016 
3017  // Allocate Healpix projection
3018  m_proj = new GHealpix;
3019 
3020  // Read WCS information from FITS header
3021  m_proj->read(table);
3022 
3023  // Set number of pixels based on NSIDE parameter
3024  m_num_pixels = static_cast<GHealpix*>(m_proj)->npix();
3025  #if defined(G_READ_HEALPIX_DEBUG)
3026  std::cout << "m_num_pixels=" << m_num_pixels << std::endl;
3027  #endif
3028 
3029  // Number of map pixels has to be a multiple of the number of
3030  // rows in column
3031  if (m_num_pixels % nrows != 0) {
3032  std::string msg = "Number of "+gammalib::str(nrows)+" rows in FITS "
3033  "table must correspond to number of "+
3034  gammalib::str(m_num_pixels)+" pixels in HealPix "
3035  "projection.";
3037  }
3038 
3039  // Determine vector length for HEALPix data storage
3040  int nentry = m_num_pixels / nrows;
3041  #if defined(G_READ_HEALPIX_DEBUG)
3042  std::cout << "nentry=" << nentry << std::endl;
3043  #endif
3044 
3045  // Determine number of maps from the number of maps that fit into
3046  // all columns. Only count columns that can fully hold the map.
3047  m_num_maps = 0;
3048  for (int icol = 0; icol < ncols; ++icol) {
3049  const GFitsTableCol* col = table[icol];
3050  if (col->number() % nentry == 0) {
3051  m_num_maps += col->number() / nentry;
3052  }
3053  }
3054  #if defined(G_READ_HEALPIX_DEBUG)
3055  std::cout << "m_num_maps=" << m_num_maps << std::endl;
3056  #endif
3057 
3058  // Initialise shape of maps
3059  m_shape.clear();
3060  m_shape.push_back(m_num_maps);
3061 
3062  // Allocate pixels to hold the map
3064 
3065  // Initialise map counter
3066  int imap = 0;
3067 
3068  // Loop over all columns
3069  for (int icol = 0; icol < ncols; ++icol) {
3070 
3071  // Get next column
3072  const GFitsTableCol* col = table[icol];
3073 
3074  // Only consider columns that can fully hold maps
3075  if (col->number() % nentry == 0) {
3076 
3077  // Determine number of maps in column
3078  int num = col->number() / nentry;
3079 
3080  // Loop over all maps in column
3081  int inx_start = 0;
3082  int inx_end = nentry;
3083  for (int i = 0; i < num; ++i) {
3084 
3085  // Load map
3086  double *ptr = m_pixels.data() + m_num_pixels*imap;
3087  for (int row = 0; row < col->nrows(); ++row) {
3088  for (int inx = inx_start; inx < inx_end; ++inx) {
3089  *ptr++ = col->real(row,inx);
3090  }
3091  }
3092  #if defined(G_READ_HEALPIX_DEBUG)
3093  std::cout << "Load map=" << imap << " index="
3094  << inx_start << "-" << inx_end << std::endl;
3095  #endif
3096 
3097  // Increment index range
3098  inx_start = inx_end;
3099  inx_end += nentry;
3100 
3101  // Increment map counter
3102  imap++;
3103 
3104  // Break if we have loaded all maps
3105  if (imap >= m_num_maps) {
3106  break;
3107  }
3108 
3109  } // endfor: looped over all maps in column
3110  } // endif: column could fully hold maps
3111 
3112  // Break if we have loaded all maps
3113  if (imap >= m_num_maps) {
3114  break;
3115  }
3116 
3117  } // endfor: looped over all columns
3118 
3119  // Return
3120  return;
3121 }
3122 
3123 
3124 /***********************************************************************//**
3125  * @brief Read WCS image from FITS HDU
3126  *
3127  * @param[in] image FITS image.
3128  *
3129  * @exception GException::invalid_argument
3130  * FITS image has less than two dimensions
3131  * @exception GException::invalid_value
3132  * Sky map covers more than 360 deg in longitude or 180 deg in
3133  * latitude
3134  *
3135  * Reads sky maps from a FITS image extension containing a set of maps given
3136  * in World Coordinate System. The method handles general n-dimensional
3137  * images and sets the map shape attribute according to the number of map
3138  * dimensions found in the FITS HDU.
3139  ***************************************************************************/
3140 void GSkyMap::read_wcs(const GFitsImage& image)
3141 {
3142  // Throw an exception if the FITS image is not at least a 2D image
3143  if (image.naxis() < 2) {
3144  std::string msg = "Sky map has "+gammalib::str(image.naxis())+
3145  " dimensions, which is less than the two dimensions"
3146  " that are required for a WCS image.";
3148  }
3149 
3150  // Allocate WCS
3151  alloc_wcs(image);
3152 
3153  // Read projection information from FITS header
3154  m_proj->read(image);
3155 
3156  // Clear map shape
3157  m_shape.clear();
3158 
3159  // Extract map dimension, number of maps and map shape from FITS image
3160  m_num_x = image.naxes(0);
3161  m_num_y = image.naxes(1);
3162  if (image.naxis() == 2) {
3163  m_num_maps = 1;
3164  m_shape.push_back(m_num_maps);
3165  }
3166  else {
3167  m_num_maps = image.naxes(2);
3168  m_shape.push_back(m_num_maps);
3169  for (int i = 3; i < image.naxis(); ++i) {
3170  m_num_maps *= image.naxes(i);
3171  m_shape.push_back(image.naxes(i));
3172  }
3173  }
3174  #if defined(G_READ_WCS_DEBUG)
3175  std::cout << "m_num_x=" << m_num_x << std::endl;
3176  std::cout << "m_num_y=" << m_num_y << std::endl;
3177  std::cout << "m_num_maps=" << m_num_maps << std::endl;
3178  #endif
3179 
3180  // Compute number of pixels
3182  #if defined(G_READ_WCS_DEBUG)
3183  std::cout << "m_num_pixels=" << m_num_pixels << std::endl;
3184  #endif
3185 
3186  // Allocate pixels to hold the map
3187  m_pixels = GNdarray(m_num_x, m_num_y, m_num_maps);
3188 
3189  // Read image
3190  double* ptr = m_pixels.data();
3191  int size = m_num_pixels * m_num_maps;
3192  for (int i = 0; i < size; ++i) {
3193  *ptr++ = (double)image.pixel(i);
3194  }
3195 
3196  // Check if the map is too large
3197  double x_size = std::abs(m_num_x * static_cast<GWcs*>(m_proj)->cdelt(0));
3198  double y_size = std::abs(m_num_y * static_cast<GWcs*>(m_proj)->cdelt(1));
3199  if (x_size > 360.001) {
3200  std::string msg = "Skymap covers "+gammalib::str(x_size)+" degrees "
3201  "in the X axis which results in ambiguous "
3202  "coordinate transformations. Please provide a "
3203  "skymap that does not cover more than 360 degrees.";
3205  }
3206  if (y_size > 180.001) {
3207  std::string msg = "Skymap covers "+gammalib::str(y_size)+" degrees "
3208  "in the Y axis which results in ambiguous "
3209  "coordinate transformations. Please provide a "
3210  "skymap that does not cover more than 180 degrees.";
3212  }
3213 
3214  // Return
3215  return;
3216 }
3217 
3218 
3219 /***********************************************************************//**
3220  * @brief Allocate WCS class
3221  *
3222  * @param[in] image FITS image.
3223  *
3224  * @exception GException::invalid_argument
3225  * CTYPE1 and CTYPE2 keywords are incompatible or WCS projection
3226  * of FITS file not supported by GammaLib.
3227  ***************************************************************************/
3228 void GSkyMap::alloc_wcs(const GFitsImage& image)
3229 {
3230  // Get standard keywords
3231  std::string ctype1 = image.string("CTYPE1");
3232  std::string ctype2 = image.string("CTYPE2");
3233 
3234  // Extract projection type
3235  std::string xproj = ctype1.substr(5,3);
3236  std::string yproj = ctype2.substr(5,3);
3237 
3238  // Check that projection type is identical on both axes
3239  if (xproj != yproj) {
3240  std::string msg = "X axis projection \""+xproj+"\" differs from Y "
3241  "axis projection \""+yproj+"\". Please specify a "
3242  "FITS image with a valid WCS header.";
3244  }
3245 
3246  // Allocate WCS registry
3247  GWcsRegistry registry;
3248 
3249  // Allocate projection from registry
3250  m_proj = registry.alloc(xproj);
3251 
3252  // Signal if projection type is not known
3253  if (m_proj == NULL) {
3254  std::string msg = "WCS projection code \""+xproj+"\" not known. Please "
3255  "specify a FITS image with one of the following "
3256  "WCS projections: "+registry.content()+".";
3258  }
3259 
3260  // Return
3261  return;
3262 }
3263 
3264 
3265 /***********************************************************************//**
3266  * @brief Create FITS HDU containing Healpix data
3267  *
3268  * This method allocates a binary table HDU that contains the Healpix data.
3269  * Deallocation of the table has to be done by the client.
3270  ***************************************************************************/
3272 {
3273  // Initialise result to NULL pointer
3274  GFitsBinTable* hdu = NULL;
3275 
3276  // Compute size of Healpix data
3277  int size = m_num_pixels * m_num_maps;
3278 
3279  // Continue only if we have pixels
3280  if (size > 0) {
3281 
3282  // Set number of rows and columns
3283  int rows = m_num_pixels;
3284  int number = m_num_maps;
3285 
3286  // Create column to hold Healpix data
3287  GFitsTableDoubleCol column = GFitsTableDoubleCol("DATA", rows, number);
3288 
3289  // Fill data into column
3290  const double* ptr = m_pixels.data();
3291  for (int inx = 0; inx < number; ++inx) {
3292  for (int row = 0; row < rows; ++row) {
3293  column(row,inx) = *ptr++;
3294  }
3295  }
3296 
3297  // Create HDU that contains Healpix map in a binary table
3298  hdu = new GFitsBinTable(rows);
3299  hdu->append(column);
3300 
3301  } // endif: there were pixels
3302 
3303  // ... otherwise create an empty header
3304  else {
3305  hdu = new GFitsBinTable;
3306  }
3307 
3308  // Set extension name
3309  hdu->extname("HEALPIX");
3310 
3311  // If we have WCS information then write into FITS header
3312  if (m_proj != NULL) {
3313  m_proj->write(*hdu);
3314  }
3315 
3316  // Set additional keywords
3317  hdu->card("NBRBINS", m_num_maps, "Number of HEALPix maps");
3318 
3319  // Return HDU
3320  return hdu;
3321 }
3322 
3323 
3324 /***********************************************************************//**
3325  * @brief Create FITS HDU containing WCS image
3326  *
3327  * This method allocates an image HDU that contains the WCS image data.
3328  * Deallocation of the image has to be done by the client.
3329  *
3330  * @todo Set additional keywords.
3331  ***************************************************************************/
3333 {
3334  // Initialise result to NULL pointer
3335  GFitsImageDouble* hdu = NULL;
3336 
3337  // Compute size of Healpix data
3338  int size = m_num_pixels * m_num_maps;
3339 
3340  // Continue only if we have pixels
3341  if (size > 0) {
3342 
3343  // Set dimension vector of all axes. In case that only one map
3344  // exists then create simply a 2D image
3345  std::vector<int> naxes;
3346  naxes.push_back(m_num_x);
3347  naxes.push_back(m_num_y);
3348  if (m_num_maps > 1) {
3349  for (int i = 0; i < ndim(); ++i) {
3350  naxes.push_back(m_shape[i]);
3351  }
3352  }
3353 
3354  // Allocate image
3355  hdu = new GFitsImageDouble(naxes);
3356 
3357  // Store data in image
3358  if (naxes.size() == 2) {
3359  const double* ptr = m_pixels.data();
3360  for (int iy = 0; iy < m_num_y; ++iy) {
3361  for (int ix = 0; ix < m_num_x; ++ix) {
3362  (*hdu)(ix,iy) = *ptr++;
3363  }
3364  }
3365  }
3366  else {
3367  const double* ptr = m_pixels.data();
3368  for (int imap = 0; imap < m_num_maps; ++imap) {
3369  for (int iy = 0; iy < m_num_y; ++iy) {
3370  for (int ix = 0; ix < m_num_x; ++ix) {
3371  (*hdu)(ix,iy,imap) = *ptr++;
3372  }
3373  }
3374  }
3375  }
3376 
3377  } // endif: there were pixels
3378 
3379  // ... otherwise create an empty header
3380  else {
3381  hdu = new GFitsImageDouble;
3382  }
3383 
3384  // Set extension name
3385  hdu->extname("IMAGE");
3386 
3387  // If we have sky projection information then write into FITS header
3388  if (m_proj != NULL) {
3389  m_proj->write(*hdu);
3390  }
3391 
3392  // Set additional keywords
3393  //TODO
3394 
3395  // Return HDU
3396  return hdu;
3397 }
3398 
3399 
3400 /***********************************************************************//**
3401  * @brief Compute solid angle subtended by 4 sky directions
3402  *
3403  * @param[in] dir1 First sky direction.
3404  * @param[in] dir2 Second sky direction.
3405  * @param[in] dir3 Third sky direction.
3406  * @param[in] dir4 Forth sky direction.
3407  * @return Solid angle (steradians).
3408  *
3409  * Estimate the solid angle subtended by 4 sky directions using Huilier's
3410  * theorem.
3411  *
3412  * Below, the definiton of the pixel cornes and sides are shown as used
3413  * within the code.
3414  *
3415  * a12
3416  * 1---------2
3417  * |\ /|
3418  * | \a13 / |
3419  * | \ / |
3420  * | \ / |
3421  * a14| X |a23
3422  * | / \ |
3423  * | / \ |
3424  * | /a24 \ |
3425  * |/ \|
3426  * 4---------3
3427  * a34
3428  *
3429  ***************************************************************************/
3430 double GSkyMap::solidangle(const GSkyDir& dir1, const GSkyDir& dir2,
3431  const GSkyDir& dir3, const GSkyDir& dir4) const
3432 {
3433  // Initialise solid angle
3434  double solidangle = 0.0;
3435 
3436  // Compute angular distances between pixel corners
3437  double a12 = dir1.dist(dir2);
3438  double a14 = dir1.dist(dir4);
3439  double a23 = dir2.dist(dir3);
3440  double a24 = dir2.dist(dir4);
3441  double a34 = dir3.dist(dir4);
3442 
3443  // Special case: a12 or a14 is zero, then pixel is a triangle composed
3444  // of [2,3,4]
3445  if (a12 <= 0.0 || a14 <= 0.0) {
3446  double s = 0.5 * (a23 + a34 + a24);
3447  solidangle = 4.0 * std::atan(std::sqrt(std::tan(0.5*s) *
3448  std::tan(0.5*(s-a23)) *
3449  std::tan(0.5*(s-a34)) *
3450  std::tan(0.5*(s-a24))));
3451  }
3452 
3453  // Special case: a23 or a 34 is zero, then pixel is a triangle composed
3454  // of [1,2,4]
3455  else if (a23 <= 0.0 || a34 <= 0.0) {
3456  double s = 0.5 * (a12 + a24 + a14);
3457  solidangle = 4.0 * std::atan(std::sqrt(std::tan(0.5*s) *
3458  std::tan(0.5*(s-a12)) *
3459  std::tan(0.5*(s-a24)) *
3460  std::tan(0.5*(s-a14))));
3461  }
3462 
3463  // Otherwise we have a polygon
3464  else {
3465 
3466  // Triangle 1 [1,2,4]
3467  double s1 = 0.5 * (a12 + a24 + a14);
3468  double excess1 = std::atan(std::sqrt(std::tan(0.5*s1) *
3469  std::tan(0.5*(s1-a12)) *
3470  std::tan(0.5*(s1-a24)) *
3471  std::tan(0.5*(s1-a14))));
3472 
3473  // Triangle 2 [2,3,4]
3474  double s2 = 0.5 * (a23 + a34 + a24);
3475  double excess2 = std::atan(std::sqrt(std::tan(0.5*s2) *
3476  std::tan(0.5*(s2-a23)) *
3477  std::tan(0.5*(s2-a34)) *
3478  std::tan(0.5*(s2-a24))));
3479 
3480  // Determine solid angle
3481  solidangle = 4.0 * (excess1 + excess2);
3482 
3483  } // endif: we had a polynom
3484 
3485  // Return solid angle
3486  return solidangle;
3487 }
3488 
3489 
3490 /***********************************************************************//**
3491  * @brief Compute solid angle subtended by 3 sky directions
3492  *
3493  * @param[in] dir1 First sky direction.
3494  * @param[in] dir2 Second sky direction.
3495  * @param[in] dir3 Third sky direction.
3496  * @return Solid angle (steradians).
3497  *
3498  * Estimate the solid angle subtended by 3 sky directions using Huilier's
3499  * theorem.
3500  *
3501  * Below, the definiton of the pixel cornes and sides are shown as used
3502  * within the code.
3503  *
3504  * a12
3505  * 1---------2
3506  * | /
3507  * | /
3508  * | /
3509  * | /
3510  * a13| /a23
3511  * | /
3512  * | /
3513  * | /
3514  * |/
3515  * 3
3516  *
3517  ***************************************************************************/
3518 double GSkyMap::solidangle(const GSkyDir& dir1, const GSkyDir& dir2,
3519  const GSkyDir& dir3) const
3520 {
3521  // Compute angular distances between pixel corners
3522  double a12 = dir1.dist(dir2);
3523  double a13 = dir1.dist(dir3);
3524  double a23 = dir2.dist(dir3);
3525 
3526  // Compute solid angle
3527  double s = 0.5 * (a12 + a23 + a13);
3528  double solidangle = 4.0 * std::atan(std::sqrt(std::tan(0.5*s) *
3529  std::tan(0.5*(s-a12)) *
3530  std::tan(0.5*(s-a23)) *
3531  std::tan(0.5*(s-a13))));
3532 
3533  // Return solid angle
3534  return solidangle;
3535 }
3536 
3537 
3538 /***********************************************************************//**
3539  * @brief Checks whether a circular region overlaps with this map
3540  *
3541  * @param[in] region Circular region
3542  * @return True if circular region overlaps with map, false otherwise
3543  *
3544  * The check is done by first testing whether the central pointing
3545  * position of the observation falls into the sky map. If this is false,
3546  * then pixel positions that border the sky map are tested for whether or
3547  * not they fall into the observation region. The positions tested can be
3548  * visualized as follows, where '*' marks the positions tested.
3549  *
3550  * * * * * * *
3551  * +---+---+---+---+
3552  * * |0,2|1,2|2,2|3,2| *
3553  * +---+---+---+---+
3554  * * |0,1|1,1|2,1|3,1| *
3555  * +---+---+---+---+
3556  * * |0,0|1,0|2,0|3,0| *
3557  * +---+---+---+---+
3558  * * * * * * *
3559  *
3560  ***************************************************************************/
3562 {
3563  // Initialise overlap with true
3564  bool overlap = false;
3565 
3566  // If the centre of the region is inside the map then signal overlap ...
3567  if (contains(region.centre())) {
3568  overlap = true;
3569  }
3570 
3571  // ... otherwise loop through each of the outer bins and check if the
3572  // region overlaps with them (within some delta)
3573  else {
3574 
3575  // Initialise test pixels
3576  GSkyPixel pix1(0.0,0.0);
3577  GSkyPixel pix2(0.0,0.0);
3578 
3579  // Check the bottom and top bins
3580  pix1.y(-1.0);
3581  pix2.y(double(ny()));
3582  for (double xbin = -1.0; xbin <= nx(); xbin += 1.0) {
3583  pix1.x(xbin);
3584  pix2.x(xbin);
3585  if ((region.contains(pix2dir(pix1))) ||
3586  (region.contains(pix2dir(pix2)))) {
3587  overlap = true;
3588  break;
3589  }
3590  } // endfor: loop on xbin
3591 
3592  // If no overlap has been found then check now the left and right
3593  // bins
3594  if (!overlap) {
3595  pix1.x(-1.0);
3596  pix2.x(double(nx()));
3597  for (double ybin = 0.0; ybin < ny(); ybin += 1.0) {
3598  pix1.y(ybin);
3599  pix2.y(ybin);
3600  if ((region.contains(pix2dir(pix1))) ||
3601  (region.contains(pix2dir(pix2)))) {
3602  overlap = true;
3603  break;
3604  }
3605  } // endfor: loop on ybin
3606  }
3607 
3608  } // endelse: check border bins
3609 
3610  // Return overlap flag
3611  return overlap;
3612 }
3613 
3614 
3615 /***********************************************************************//**
3616  * @brief Check if HDU contains HEALPix data
3617  *
3618  * @param[in] hdu FITS Header Data Unit.
3619  * @return True is HDU contains HEALPix data.
3620  *
3621  * Returns true if the HDU is not an image and the HDU has the "PIXTYPE"
3622  * keyword set to "HEALPIX".
3623  ***************************************************************************/
3624 bool GSkyMap::is_healpix(const GFitsHDU& hdu) const
3625 {
3626  // Initialise flag
3627  bool flag(false);
3628 
3629  // If PIXTYPE keyword equals "HEALPIX" then signal that we have
3630  // HEALPix data
3631  if ((hdu.exttype() != GFitsHDU::HT_IMAGE) &&
3632  (hdu.has_card("PIXTYPE") && (hdu.string("PIXTYPE") == "HEALPIX"))) {
3633  flag = true;
3634  }
3635 
3636  // Return flag
3637  return (flag);
3638 }
3639 
3640 
3641 /***********************************************************************//**
3642  * @brief Check if HDU contains WCS data
3643  *
3644  * @param[in] hdu FITS Header Data Unit.
3645  * @return True is HDU contains WCS data.
3646  *
3647  * Returns true if the HDU is an image and the HDU has the "NAXIS" keyword
3648  * set to a value >= 2.
3649  ***************************************************************************/
3650 bool GSkyMap::is_wcs(const GFitsHDU& hdu) const
3651 {
3652  // Initialise flag
3653  bool flag(false);
3654 
3655  // If extension is an image thn signal that we have WCS data
3656  if ((hdu.exttype() == GFitsHDU::HT_IMAGE) &&
3657  (hdu.has_card("NAXIS") && (hdu.integer("NAXIS") >= 2))) {
3658  flag = true;
3659  }
3660 
3661  // Return flag
3662  return (flag);
3663 }
3664 
3665 
3666 /***********************************************************************//**
3667  * @brief Check if map is the same
3668  *
3669  * @param[in] map Sky map.
3670  * @return True is sky map definition is identical.
3671  *
3672  * Returns true if the sky map definition is identical to the one of the
3673  * current map.
3674  ***************************************************************************/
3675 bool GSkyMap::is_same(const GSkyMap& map) const
3676 {
3677  // Initialise flag
3678  bool flag = true;
3679 
3680  // Single loop that can be exited any time
3681  do {
3682 
3683  // If sky projection differs then set flag to false and break
3684  if (m_proj->code() != map.m_proj->code()) {
3685  flag = false;
3686  break;
3687  }
3688 
3689  // If sky projection has different coordinate system then set flag
3690  // to false and break
3691  if (m_proj->coordsys() != map.m_proj->coordsys()) {
3692  flag = false;
3693  break;
3694  }
3695 
3696  // If sky map dimension differs then set flag to false and break
3697  if ((m_num_x != map.m_num_x) || (m_num_y != map.m_num_y)) {
3698  flag = false;
3699  break;
3700  }
3701 
3702  // For WCS projections, if sky maps have different reference values,
3703  // pixels or pixel sizes then set flag to false and break
3704  const GWcs* wcs1 = dynamic_cast<const GWcs*>(m_proj);
3705  const GWcs* wcs2 = dynamic_cast<const GWcs*>(map.m_proj);
3706  if ((wcs1 != NULL) && (wcs2 != NULL)) {
3707  if ((wcs1->crval(0) != wcs2->crval(0)) ||
3708  (wcs1->crval(1) != wcs2->crval(1)) ||
3709  (wcs1->crpix(0) != wcs2->crpix(0)) ||
3710  (wcs1->crpix(1) != wcs2->crpix(1)) ||
3711  (wcs1->cdelt(0) != wcs2->cdelt(0)) ||
3712  (wcs1->cdelt(1) != wcs2->cdelt(1))) {
3713  flag = false;
3714  break;
3715  }
3716  }
3717 
3718  // For Healpix projections, if sky maps have different number of pixels,
3719  // nside or ordering then set flag to false and break
3720  const GHealpix* healpix1 = dynamic_cast<const GHealpix*>(m_proj);
3721  const GHealpix* healpix2 = dynamic_cast<const GHealpix*>(map.m_proj);
3722  if ((healpix1 != NULL) && (healpix2 != NULL)) {
3723  if ((healpix1->npix() != healpix2->npix()) ||
3724  (healpix1->nside() != healpix2->nside()) ||
3725  (healpix1->ordering() != healpix2->ordering())) {
3726  flag = false;
3727  break;
3728  }
3729  }
3730 
3731  } while(false);
3732 
3733  // Return flag
3734  return (flag);
3735 }
3736 
3737 
3738 /***********************************************************************//**
3739  * @brief Convolve sky map
3740  *
3741  * @param[in] kernel Convolution kernel type ("DISK", "GAUSSIAN").
3742  * @param[in] par Convolution parameter.
3743  * @param[in] normalise Normalise kernel?
3744  *
3745  * Convolves all sky maps using the specified @p kernel and a convolution
3746  * parameter.
3747  *
3748  * For the "DISK" kernel the convolution parameter is the disk radius in
3749  * degrees. For the "GAUSSIAN" kernel the convolution parameter is the
3750  * Gaussian sigma in degrees.
3751  *
3752  * If @p normalise is true the kernel is normalised to an integral of unity
3753  * so that the operation corresponds to a smoothing. Otherwise the maximum
3754  * of the kernel is 1 so that the convolution is a correlation.
3755  ***************************************************************************/
3756 void GSkyMap::convolve(const std::string& kernel,
3757  const double& par,
3758  const bool& normalise)
3759 {
3760  // Continue only if the convolution parameter is positive
3761  if (par > 0.0) {
3762 
3763  // Get FFT of convolution kernel
3764  GFft fft_kernel(convolution_kernel(kernel, par, normalise));
3765 
3766  // Loop over all sky maps
3767  for (int k = 0; k < m_num_maps; ++k) {
3768 
3769  // Extract sky map
3770  GNdarray array(m_num_x, m_num_y);
3771  const double *src = m_pixels.data() + k*m_num_pixels;
3772  double *dst = array.data();
3773  for (int i = 0; i < m_num_pixels; ++i) {
3774  *dst++ = *src++;
3775  }
3776 
3777  // FFT of sky map
3778  GFft fft_array = GFft(array);
3779 
3780  // Multiply FFT of sky map with FFT of kernel
3781  GFft fft_smooth = fft_array * fft_kernel;
3782 
3783  // Backward transform sky map
3784  GNdarray smooth = fft_smooth.backward();
3785 
3786  // Insert sky map
3787  src = smooth.data();
3788  dst = m_pixels.data() + k*m_num_pixels;
3789  for (int i = 0; i < m_num_pixels; ++i) {
3790  *dst++ = *src++;
3791  }
3792 
3793  } // endfor: looped over all sky map
3794 
3795  } // endif: convolution parameter was positive
3796 
3797  // Return
3798  return;
3799 }
3800 
3801 /***********************************************************************//**
3802  * @brief Return convolution kernel
3803  *
3804  * @param[in] kernel Convolution kernel type ("DISK", "GAUSSIAN").
3805  * @param[in] par Convolution parameter (>0).
3806  * @param[in] normalise Normalise kernel?
3807  * @return Array filled with convolution kernel.
3808  *
3809  * If @p normalise is true the kernel is normalised to an integral of unity
3810  * so that the operation corresponds to a smoothing. Otherwise the maximum
3811  * of the kernel is 1 so that the convolution is a correlation.
3812  ***************************************************************************/
3813 GNdarray GSkyMap::convolution_kernel(const std::string& kernel,
3814  const double& par,
3815  const bool& normalise) const
3816 {
3817  // Get pointer on WCS projection
3818  const GWcs* wcs = dynamic_cast<const GWcs*>(m_proj);
3819  if (wcs == NULL) {
3820  std::string msg = "Sky map is not a WCS projection. Method is only "
3821  "valid for WCS projections.";
3823  }
3824 
3825  // Get X and Y step size
3826  double dx = wcs->cdelt(0);
3827  double dy = wcs->cdelt(1);
3828 
3829  // Initialise kernel
3830  GNdarray kern(m_num_x, m_num_y);
3831 
3832  // Initialise kernel sum
3833  double sum = 0.0;
3834 
3835  // Handle disk kernel
3836  if (gammalib::toupper(kernel) == "DISK") {
3837 
3838  // Fill kernel
3839  for (int ix1 = 0, ix2 = m_num_x; ix1 < m_num_x; ++ix1, --ix2) {
3840  double x = ix1 * dx;
3841  double xqs = x * x;
3842  for (int iy1 = 0, iy2 = m_num_y; iy1 < m_num_y; ++iy1, --iy2) {
3843  double y = iy1 * dy;
3844  if (std::sqrt(xqs + y*y) <= par) {
3845  kern(ix1,iy1) += 1.0;
3846  sum += 1.0;
3847  if (ix2 < m_num_x) {
3848  kern(ix2,iy1) += 1.0;
3849  sum += 1.0;
3850  }
3851  if (iy2 < m_num_y) {
3852  kern(ix1,iy2) += 1.0;
3853  sum += 1.0;
3854  }
3855  if ((ix2 < m_num_x) && (iy2 < m_num_y)) {
3856  kern(ix2,iy2) += 1.0;
3857  sum += 1.0;
3858  }
3859  }
3860  }
3861  }
3862 
3863  } // endif: handled disk kernel
3864 
3865  // Handle Gaussian kernel
3866  else if (gammalib::toupper(kernel) == "GAUSSIAN") {
3867 
3868  // Compute the exponential normalisation
3869  double norm = -0.5 / (par * par);
3870 
3871  // Fill kernel
3872  for (int ix1 = 0, ix2 = m_num_x; ix1 < m_num_x; ++ix1, --ix2) {
3873  double x = ix1 * dx;
3874  double xqs = x * x;
3875  for (int iy1 = 0, iy2 = m_num_y; iy1 < m_num_y; ++iy1, --iy2) {
3876  double y = iy1 * dy;
3877  double value = std::exp(norm * (xqs + y*y));
3878  kern(ix1,iy1) += value;
3879  sum += value;
3880  if (ix2 < m_num_x) {
3881  kern(ix2,iy1) += value;
3882  sum += value;
3883  }
3884  if (iy2 < m_num_y) {
3885  kern(ix1,iy2) += value;
3886  sum += value;
3887  }
3888  if ((ix2 < m_num_x) && (iy2 < m_num_y)) {
3889  kern(ix2,iy2) += value;
3890  sum += value;
3891  }
3892  }
3893  }
3894 
3895  } // endif: handled Gaussian kernel
3896 
3897  // ... otherwise throw an exception
3898  else {
3899  std::string msg = "Invalid kernel type \""+kernel+"\" specified. "
3900  "Please specify one of \"DISK\", \"GAUSSIAN\".";
3902  }
3903 
3904  // Normalise kernel
3905  if ((normalise) && (sum > 0.0)) {
3906  for (int ix = 0; ix < m_num_x; ++ix) {
3907  for (int iy = 0; iy < m_num_y; ++iy) {
3908  kern(ix,iy) /= sum;
3909  }
3910  }
3911  }
3912 
3913  // Return kernel
3914  return kern;
3915 }
3916 
3917 
3918 /*==========================================================================
3919  = =
3920  = Friends =
3921  = =
3922  ==========================================================================*/
3923 
3924 /***********************************************************************//**
3925  * @brief Computes square root of sky map elements
3926  *
3927  * @param[in] map Sky map.
3928  * @return Sky map containing the square root of every element.
3929  ***************************************************************************/
3930 GSkyMap sqrt(const GSkyMap& map)
3931 {
3932  // Initialise result vector
3933  GSkyMap result(map);
3934 
3935  // Loop over all maps
3936  for (int i = 0; i < map.nmaps(); ++i) {
3937 
3938  // Loop over all bins
3939  for (int j = 0; j < map.npix(); ++j) {
3940 
3941  // Get the content from the bin
3942  double content = map(j,i);
3943 
3944  // Throw an exception if content is negative
3945  if (content < 0.0) {
3946  std::string msg = "Negative value encountered. "
3947  "Cannot take the sqrt from a negative "
3948  "value";
3949  throw GException::invalid_value(G_SQRT, msg);
3950  }
3951 
3952  // Set content of the result map
3953  result(j,i) = std::sqrt(content);
3954 
3955  } // endfor: looped over all bins
3956 
3957  } // endfor: looped over all maps
3958 
3959  // Return sky map
3960  return result;
3961 }
3962 
3963 
3964 /***********************************************************************//**
3965  * @brief Computes the natural logarithm of sky map elements
3966  *
3967  * @param[in] map Sky map.
3968  * @return Sky map containing the natural logarithm of every element.
3969  ***************************************************************************/
3970 GSkyMap log(const GSkyMap& map)
3971 {
3972  // Initialise result vector
3973  GSkyMap result(map);
3974 
3975  // Loop over all maps
3976  for (int i = 0; i < map.nmaps(); ++i) {
3977 
3978  // Loop over all bins
3979  for (int j = 0; j < map.npix(); ++j) {
3980 
3981  // Get the content from the bin
3982  double content = map(j,i);
3983 
3984  // Throw an exception if content is negative
3985  if (content <= 0.0) {
3986  std::string msg = "Negative or zero value encountered. "
3987  "Cannot calculate the logarithm of a "
3988  "negative or zero value";
3989  throw GException::invalid_value(G_LOG, msg);
3990  }
3991 
3992  // Set content of the result map
3993  result(j,i) = std::log(content);
3994 
3995  } // endfor: looped over all bins
3996 
3997  } // endfor: looped over all maps
3998 
3999  // Return sky map
4000  return result;
4001 }
4002 
4003 
4004 /***********************************************************************//**
4005  * @brief Computes the base 10 logarithm of sky map elements
4006  *
4007  * @param[in] map Sky map.
4008  * @return Sky map containing the base 10 logarithm of every element.
4009  ***************************************************************************/
4010 GSkyMap log10(const GSkyMap& map)
4011 {
4012  // Initialise result vector
4013  GSkyMap result(map);
4014 
4015  // Loop over all maps
4016  for (int i = 0; i < map.nmaps(); ++i) {
4017 
4018  // Loop over all bins
4019  for (int j = 0; j < map.npix(); ++j) {
4020 
4021  // Get the content from the bin
4022  double content = map(j,i);
4023 
4024  // Throw an exception if content is negative
4025  if (content <= 0.0) {
4026  std::string msg = "Negative or zero value encountered. "
4027  "Cannot calculate the logarithm of a "
4028  "negative or zero value";
4029  throw GException::invalid_value(G_LOG10, msg);
4030  }
4031 
4032  // Set content of the result map
4033  result(j,i) = std::log10(content);
4034 
4035  } // endfor: looped over all bins
4036 
4037  } // endfor: looped over all maps
4038 
4039  // Return sky map
4040  return result;
4041 }
4042 
4043 
4044 /***********************************************************************//**
4045  * @brief Computes the absolute value of sky map elements
4046  *
4047  * @param[in] map Sky map.
4048  * @return Sky map containing the absolute value of every element.
4049  ***************************************************************************/
4050 GSkyMap abs(const GSkyMap& map)
4051 {
4052  // Initialise result vector
4053  GSkyMap result(map);
4054 
4055  // Loop over all maps
4056  for (int i = 0; i < map.nmaps(); ++i) {
4057 
4058  // Loop over all bins
4059  for (int j = 0; j < map.npix(); ++j) {
4060 
4061  // Get the content from the bin
4062  double content = map(j,i);
4063 
4064  // Set content of the result map
4065  result(j,i) = std::abs(content);
4066 
4067  } // endfor: looped over all bins
4068 
4069  } // endfor: looped over all maps
4070 
4071  // Return sky map
4072  return result;
4073 }
4074 
4075 
4076 /***********************************************************************//**
4077  * @brief Computes the sign value of sky map elements
4078  *
4079  * @param[in] map Sky map.
4080  * @return Sky map containing the sign value of every pixel.
4081  *
4082  * This method returns a sky map filled with a value of 1 if the pixel
4083  * is positive, a value of -1 if the pixel is negative or a value of 0
4084  * if the pixel is 0.
4085  ***************************************************************************/
4086 GSkyMap sign(const GSkyMap& map)
4087 {
4088  // Initialise result vector
4089  GSkyMap result(map);
4090 
4091  // Loop over all maps
4092  for (int i = 0; i < map.nmaps(); ++i) {
4093 
4094  // Loop over all bins
4095  for (int j = 0; j < map.npix(); ++j) {
4096 
4097  // Get the content from the bin
4098  double content = map(j,i);
4099 
4100  // Handle the 3 cases
4101  if (content < 0.0) {
4102  result(j,i) = -1.0;
4103  }
4104  else if (content > 0.0) {
4105  result(j,i) = +1.0;
4106  }
4107  else {
4108  result(j,i) = 0.0;
4109  }
4110 
4111  } // endfor: looped over all bins
4112 
4113  } // endfor: looed over all maps
4114 
4115  // Return sky map
4116  return result;
4117 }
4118 
4119 
4120 /***********************************************************************//**
4121  * @brief Clips map at given value
4122  *
4123  * @param[in] map Sky map.
4124  * @param[in] thresh Threshold value.
4125  * @return Sky map containing the value in the original map if >= thresh,
4126  * otherwise thresh
4127  ***************************************************************************/
4128 GSkyMap clip(const GSkyMap& map, const double& thresh)
4129 {
4130  // Initialise result vector
4131  GSkyMap result(map);
4132 
4133  // Loop over all maps
4134  for (int i = 0; i < map.nmaps(); ++i) {
4135 
4136  // Loop over all bins
4137  for (int j = 0; j < map.npix(); ++j) {
4138 
4139  // Get the content from the bin
4140  double content = map(j,i);
4141 
4142  // If content is below threshold, set to thresh
4143  if (content < thresh) {
4144  result(j,i) = thresh;
4145  }
4146 
4147  // ... otherwise keep content
4148  else {
4149  result(j,i) = content;
4150  }
4151 
4152  } // endfor: looped over all bins
4153 
4154  } // endfor: looped over all maps
4155 
4156  // Return sky map
4157  return result;
4158 }
4159 
4160 
4161 /*==========================================================================
4162  = =
4163  = Friends =
4164  = =
4165  ==========================================================================*/
4166 
4167 /***********************************************************************//**
4168  * @brief Equality operator
4169  *
4170  * @param[in] a First sky map.
4171  * @param[in] b Second sky map.
4172  * @return True if @p a and @p b are identical.
4173  *
4174  * Two sky maps are considered identical if they have the same projections,
4175  * coordinate definiton and number of pixels. The actual content of the map
4176  * does not need to be identical.
4177  ***************************************************************************/
4178 bool operator==(const GSkyMap &a, const GSkyMap &b)
4179 {
4180  // Return result
4181  return a.is_same(b);
4182 }
4183 
4184 
4185 /***********************************************************************//**
4186  * @brief Non-equality operator
4187  *
4188  * @param[in] a First sky map.
4189  * @param[in] b Second sky map.
4190  * @return True if @p a and @p b are not identical.
4191  *
4192  * Two sky maps are considered different if they either differ in projections,
4193  * coordinate definiton or number of pixels. The actual content of the map
4194  * does is not relevant.
4195  ***************************************************************************/
4196 bool operator!=(const GSkyMap &a, const GSkyMap &b)
4197 {
4198  // Return result
4199  return !(a == b);
4200 }
GSkyDirs boundaries(const GSkyPixel &pixel, const int &step=1) const
Return pixel boundaries.
Definition: GHealpix.cpp:711
GFitsHDU * at(const int &extno)
Get pointer to HDU.
Definition: GFits.cpp:233
bool is_same(const GSkyMap &map) const
Check if map is the same.
Definition: GSkyMap.cpp:3675
GSkyMap(void)
Void constructor.
Definition: GSkyMap.cpp:116
const int & ncols(void) const
Return number of columns in table.
Definition: GFitsTable.hpp:134
void stack_maps(void)
Stack all maps into a single map.
Definition: GSkyMap.cpp:2341
std::string print(const GChatter &chatter=NORMAL) const
Print pixel.
Definition: GSkyPixel.cpp:292
GSkyMap & operator=(const GSkyMap &map)
Assignment operator.
Definition: GSkyMap.cpp:337
Interface definition for the WCS registry class.
Sky map class.
Definition: GSkyMap.hpp:89
double dist_deg(const GSkyDir &dir) const
Compute angular distance between sky directions in degrees.
Definition: GSkyDir.hpp:286
VO client class.
Definition: GVOClient.hpp:59
Double precision FITS image class definition.
FITS table double column class interface definition.
Abstract FITS image base class.
Definition: GFitsImage.hpp:43
void number(const int &number)
Set number of elements in column.
int size(void) const
Return pixel dimension.
Definition: GSkyPixel.hpp:145
#define G_EXTRACT
Definition: GSkyMap.cpp:79
const int & naxis(void) const
Return dimension of image.
Definition: GFitsImage.hpp:156
GSkyMap & operator*=(const GSkyMap &map)
Multiplication operator.
Definition: GSkyMap.cpp:578
std::string number(const std::string &noun, const int &number)
Convert singular noun into number noun.
Definition: GTools.cpp:1167
bool has_card(const int &cardno) const
Check existence of header card.
Definition: GFitsHDU.hpp:233
const int & ny(void) const
Returns number of pixels in y coordinate.
Definition: GSkyMap.hpp:379
bool is_1D(void) const
Check if pixel is 1D.
Definition: GSkyPixel.hpp:157
GSkyPixel inx2pix(const int &index) const
Converts pixel index into sky map pixel.
Definition: GSkyMap.cpp:1301
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:864
int pix2inx(const GSkyPixel &pixel) const
Converts sky map pixel into pixel index.
Definition: GSkyMap.cpp:1407
const GSkyProjection * projection(void) const
Returns pointer to sky projection.
Definition: GSkyMap.hpp:465
const int & nside(void) const
Returns number of divisions of the side of each base pixel.
Definition: GHealpix.hpp:207
GNdarray counts(void) const
Returns array with total number of counts for count maps.
Definition: GSkyMap.cpp:1496
GSkyMap * clone(void) const
Clone sky map.
Definition: GSkyMap.cpp:1118
#define G_LOG10
Definition: GSkyMap.cpp:90
void clear(void)
Clear array.
Definition: GNdarray.cpp:527
void x(const double &x)
Set x value of sky map pixel.
Definition: GSkyPixel.hpp:204
GSkyMap extract(const int &map, const int &nmaps=1) const
Extract maps into a new sky map object.
Definition: GSkyMap.cpp:2139
#define G_OP_ACCESS_1D
Definition: GSkyMap.cpp:65
int m_num_maps
Number of maps (used for pixel allocation)
Definition: GSkyMap.hpp:230
int m_num_pixels
Number of pixels (used for pixel allocation)
Definition: GSkyMap.hpp:229
void read_healpix(const GFitsTable &table)
Read Healpix data from FITS table.
Definition: GSkyMap.cpp:3008
void publish(const GFitsHDU &hdu)
Publish FITS HDU.
Definition: GVOClient.cpp:372
bool has_extno(void) const
Signal if filename has an extension number.
Definition: GFilename.hpp:267
Sky directions container class definition.
void set(const std::string &coords, const double &crval1, const double &crval2, const double &crpix1, const double &crpix2, const double &cdelt1, const double &cdelt2)
Set World Coordinate System parameters.
Definition: GWcs.cpp:859
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1253
std::vector< int > m_shape
Shape of the maps.
Definition: GSkyMap.hpp:233
Abstract FITS extension base class.
Definition: GFitsHDU.hpp:51
#define G_OP_VALUE
Definition: GSkyMap.cpp:67
bool overlaps(const GSkyRegion &region) const
Checks whether a region overlaps with this map.
Definition: GSkyMap.cpp:2102
int m_num_x
Number of pixels in x direction (only 2D)
Definition: GSkyMap.hpp:231
#define G_FLUX3
Definition: GSkyMap.cpp:74
int & index2(void)
Access index 2.
Definition: GBilinear.hpp:114
Generic matrix class definition.
Sky regions container class definition.
const std::string & type(void) const
Return region type.
Definition: GSkyRegion.hpp:137
#define G_OP_UNARY_DIV2
Definition: GSkyMap.cpp:64
const double onethird
Definition: GMath.hpp:50
#define G_PIX2DIR
Definition: GSkyMap.cpp:69
void nrows(const int &nrows)
Set number of rows in column.
bool has_extname(void) const
Signal if filename has an extension name.
Definition: GFilename.hpp:255
virtual HDUType exttype(void) const =0
GSkyProjection * m_proj
Pointer to sky projection.
Definition: GSkyMap.hpp:234
double sum(const GVector &vector)
Computes vector sum.
Definition: GVector.cpp:944
std::string extname(const std::string &defaultname="") const
Return extension name.
Definition: GFilename.cpp:385
GFitsTableCol * append(const GFitsTableCol &column)
Append column to the table.
Definition: GFitsTable.hpp:147
double crval(const int &inx) const
Return value of reference pixel.
Definition: GWcs.cpp:887
int m_num_y
Number of pixels in y direction (only 2D)
Definition: GSkyMap.hpp:232
#define G_SQRT
Definition: GSkyMap.cpp:88
double & weight2(void)
Access weight 2.
Definition: GBilinear.hpp:162
Gammalib tools definition.
virtual void read(const GFitsHDU &hdu)=0
GSkyMap & operator+=(const GSkyMap &map)
Map addition operator.
Definition: GSkyMap.cpp:398
void y(const double &y)
Set y value of sky map pixel.
Definition: GSkyPixel.hpp:221
FITS file class.
Definition: GFits.hpp:63
void 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)
Set World Coordinate System.
Definition: GSkyMap.cpp:2946
#define G_INX2DIR
Definition: GSkyMap.cpp:68
GNdarray convolution_kernel(const std::string &kernel, const double &par, const bool &normalise) const
Return convolution kernel.
Definition: GSkyMap.cpp:3813
FITS file class interface definition.
Interface for the circular sky region class.
virtual int size(void) const =0
GNdarray backward(void) const
Backward Fast Fourier Transform.
Definition: GFft.cpp:401
virtual std::string name(void) const =0
#define G_FLUX4
Definition: GSkyMap.cpp:75
void read(const GFitsHDU &hdu)
Read skymap from FITS HDU.
Definition: GSkyMap.cpp:2664
HealPix projection class definition.
Sky map class definition.
GSkyMap & operator-=(const GSkyMap &map)
Map subtraction operator.
Definition: GSkyMap.cpp:488
const std::vector< int > & shape(void) const
Returns shape of maps.
Definition: GSkyMap.hpp:403
#define G_OVERLAPS
Definition: GSkyMap.cpp:78
std::string centre(const std::string &s, const int &n, const char &c= ' ')
Centre string to achieve a length of n characters.
Definition: GTools.cpp:1118
Abstract FITS image base class definition.
void publish(const std::string &name="") const
Publish sky map.
Definition: GSkyMap.cpp:2750
World Coordinate Projection registry class interface definition.
virtual std::string code(void) const =0
VO client class interface definition.
void index(const double &index)
Set sky map pixel index.
Definition: GSkyPixel.hpp:183
bool is_healpix(const GFitsHDU &hdu) const
Check if HDU contains HEALPix data.
Definition: GSkyMap.cpp:3624
void saveto(const GFilename &filename, const bool &clobber=false)
Saves to specified FITS file.
Definition: GFits.cpp:1293
std::string content(void) const
Return list of names in registry.
Definition: GRegistry.cpp:54
std::string print(const GChatter &chatter=NORMAL) const
Print sky map.
Definition: GSkyMap.cpp:2787
void clear(void)
Clear node array.
Definition: GBilinear.cpp:156
Sky map pixel class.
Definition: GSkyPixel.hpp:74
void copy_members(const GSkyMap &map)
Copy class members.
Definition: GSkyMap.cpp:2884
#define G_OP_UNARY_ADD
Definition: GSkyMap.cpp:60
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1358
GSkyDir pix2dir(const GSkyPixel &pixel) const
Returns sky direction of pixel.
Definition: GSkyMap.cpp:1360
double & operator()(const int &index, const int &map=0)
Pixel index access operator.
Definition: GSkyMap.cpp:789
double crpix(const int &inx) const
Return reference pixel.
Definition: GWcs.cpp:911
virtual GSkyDir pix2dir(const GSkyPixel &pixel) const
Returns sky direction of pixel.
Definition: GHealpix.cpp:384
HealPix projection class interface definition.
Definition: GHealpix.hpp:50
Abstract interface for the sky region class.
Definition: GSkyRegion.hpp:57
#define G_EXTRACT_REG
Definition: GSkyMap.cpp:81
const int & nmaps(void) const
Returns number of maps.
Definition: GSkyMap.hpp:391
GSkyDir inx2dir(const int &index) const
Returns sky direction of pixel.
Definition: GSkyMap.cpp:1331
GVector tan(const GVector &vector)
Computes tangens of vector elements.
Definition: GVector.cpp:1379
virtual bool contains(const GSkyDir &dir) const =0
Double precision FITS image class.
GFitsImageDouble * create_wcs_hdu(void) const
Create FITS HDU containing WCS image.
Definition: GSkyMap.cpp:3332
GFitsHDU * write(GFits &file, const std::string &extname="") const
Write sky map into FITS file.
Definition: GSkyMap.cpp:2697
Filename class.
Definition: GFilename.hpp:62
GWcs * alloc(const std::string &code) const
Allocate World Coordinate System of given code.
Abstract interface for FITS table column.
GNdarray m_pixels
Skymap pixels.
Definition: GSkyMap.hpp:235
Fast Fourier Transformation class.
Definition: GFft.hpp:57
virtual double pixel(const int &ix) const =0
double cdelt(const int &inx) const
Return pixel size.
Definition: GWcs.cpp:935
int & index1(void)
Access index 1.
Definition: GBilinear.hpp:102
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1274
bool m_contained
Sky direction is contained in map.
Definition: GSkyMap.hpp:239
GNdarray sign(const GNdarray &array)
Computes sign of array elements.
Definition: GNdarray.cpp:1269
int & index4(void)
Access index 4.
Definition: GBilinear.hpp:138
bool contains(const GSkyDir &dir) const
Checks if sky direction falls in map.
Definition: GSkyMap.cpp:2023
const double * data(void) const
Data access method (const version)
Definition: GNdarray.hpp:389
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
int size(void) const
Return number of elements in array.
Definition: GNdarray.hpp:294
GChatter
Definition: GTypemaps.hpp:33
void alloc_wcs(const GFitsImage &image)
Allocate WCS class.
Definition: GSkyMap.cpp:3228
#define G_CONVOLUTION_KERNEL
Definition: GSkyMap.cpp:91
int integer(const std::string &keyname) const
Return card value as integer.
Definition: GFitsHDU.hpp:436
virtual void write(GFitsHDU &hdu) const =0
#define G_LOG
Definition: GSkyMap.cpp:89
void clear(void)
Clear instance.
Definition: GSkyMap.cpp:1100
const int & npix(void) const
Returns number of pixels.
Definition: GHealpix.hpp:196
void convolve(const std::string &kernel, const double &par, const bool &normalise)
Convolve sky map.
Definition: GSkyMap.cpp:3756
Vector class interface definition.
N-dimensional array class.
Definition: GNdarray.hpp:44
const std::string & extname(void) const
Return extension name.
Definition: GFitsHDU.hpp:162
#define G_SOLIDANGLE1
Definition: GSkyMap.cpp:76
double & weight3(void)
Access weight 3.
Definition: GBilinear.hpp:174
virtual std::string coordsys(void) const
Returns coordinate system.
#define G_FLUX2
Definition: GSkyMap.cpp:73
int & index3(void)
Access index 3.
Definition: GBilinear.hpp:126
virtual std::string print(const GChatter &chatter=NORMAL) const =0
Print content of object.
const int & nrows(void) const
Return number of rows in table.
Definition: GFitsTable.hpp:119
Sky region container class.
Definition: GSkyRegions.hpp:56
#define G_OP_UNARY_DIV
Definition: GSkyMap.cpp:63
int naxes(const int &axis) const
Return dimension of an image axis.
Definition: GFitsImage.cpp:344
Fast Fourier transformation class interface definition.
#define G_EXTRACT_INT
Definition: GSkyMap.cpp:80
GSkyRegionCircle region_circle(void) const
Return sky region circle that encloses the sky map.
Definition: GSkyMap.cpp:2401
void free_members(void)
Delete class members.
Definition: GSkyMap.cpp:2911
#define G_DIR2INX
Definition: GSkyMap.cpp:70
bool m_hascache
Cache is valid.
Definition: GSkyMap.hpp:238
Abstract world coordinate system base class.
Definition: GWcs.hpp:51
double & weight4(void)
Access weight 4.
Definition: GBilinear.hpp:186
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition: GFilename.hpp:189
const GSkyDir & centre(void) const
Return circular region centre.
void clear(void)
Clear sky direction.
Definition: GSkyDir.cpp:164
GSkyDir m_last_dir
Last sky direction.
Definition: GSkyMap.hpp:240
GSkyMap & operator/=(const GSkyMap &map)
Division operator.
Definition: GSkyMap.cpp:672
GSkyPixel dir2pix(const GSkyDir &dir) const
Returns sky map pixel for a given sky direction.
Definition: GSkyMap.cpp:1472
GBilinear m_interpol
Bilinear interpolator.
Definition: GSkyMap.hpp:241
GFitsBinTable * create_healpix_hdu(void) const
Create FITS HDU containing Healpix data.
Definition: GSkyMap.cpp:3271
virtual double real(const int &row, const int &inx=0) const =0
int ndim(void) const
Returns dimension of maps.
Definition: GSkyMap.hpp:415
virtual double solidangle(const GSkyPixel &pixel) const =0
double & weight1(void)
Access weight 1.
Definition: GBilinear.hpp:150
bool is_wcs(const GFitsHDU &hdu) const
Check if HDU contains WCS data.
Definition: GSkyMap.cpp:3650
void read_wcs(const GFitsImage &image)
Read WCS image from FITS HDU.
Definition: GSkyMap.cpp:3140
FITS binary table class.
bool overlaps_circle(const GSkyRegionCircle &region) const
Checks whether a circular region overlaps with this map.
Definition: GSkyMap.cpp:3561
std::string ordering(void) const
Returns ordering parameter.
Definition: GHealpix.cpp:511
int size(void) const
Return number of HDUs in FITS file.
Definition: GFits.hpp:237
GChatter reduce(const GChatter &chatter)
Reduce chattiness by one level.
Definition: GTypemaps.hpp:65
GSkyMap clip(const GSkyMap &map, const double &thresh)
Clips map at given value.
Definition: GSkyMap.cpp:4128
#define G_NMAPS
Definition: GSkyMap.cpp:58
bool contains(const GSkyDir &dir) const
Checks if sky direction lies within region.
int extno(const int &defaultno=-1) const
Return extension number.
Definition: GFilename.cpp:410
std::string string(const std::string &keyname) const
Return card value as string.
Definition: GFitsHDU.hpp:410
Exception handler interface definition.
GFitsHDU * append(const GFitsHDU &hdu)
Append HDU to FITS file.
Definition: GFits.cpp:678
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
Definition: GSkyDir.hpp:271
std::string toupper(const std::string &s)
Convert string to upper case.
Definition: GTools.cpp:941
#define G_CONSTRUCT_HPX
Definition: GSkyMap.cpp:54
#define G_FLUX1
Definition: GSkyMap.cpp:72
FITS binary table class definition.
const int & nx(void) const
Returns number of pixels in x coordinate.
Definition: GSkyMap.hpp:363
#define G_ALLOC_WCS
Definition: GSkyMap.cpp:87
#define G_OP_UNARY_SUB
Definition: GSkyMap.cpp:61
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1232
void save(const GFilename &filename, const bool &clobber=false) const
Save sky map into FITS file.
Definition: GSkyMap.cpp:2611
void smooth(const std::string &kernel, const double &par)
Smooth sky map.
Definition: GSkyMap.hpp:433
double solidangle(const int &index) const
Returns solid angle of pixel.
Definition: GSkyMap.cpp:1858
bool operator==(const GEnergy &a, const GEnergy &b)
Energy equality operator friend.
Definition: GEnergy.hpp:297
Abstract sky projection base class.
virtual GSkyDir pix2dir(const GSkyPixel &pixel) const =0
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
virtual ~GSkyMap(void)
Destructor.
Definition: GSkyMap.cpp:313
#define G_DIR2PIX
Definition: GSkyMap.cpp:71
bool is_2D(void) const
Check if pixel is 2D.
Definition: GSkyPixel.hpp:169
#define G_SOLIDANGLE2
Definition: GSkyMap.cpp:77
#define G_SET_WCS
Definition: GSkyMap.cpp:83
void load(const GFilename &filename)
Load skymap from FITS file.
Definition: GSkyMap.cpp:2511
const int & npix(void) const
Returns number of pixels.
Definition: GSkyMap.hpp:347
#define G_OP_UNARY_MUL
Definition: GSkyMap.cpp:62
GFitsHeaderCard & card(const int &cardno)
Return header card.
Definition: GFitsHDU.hpp:259
void init_members(void)
Initialise class members.
Definition: GSkyMap.cpp:2857
Sky direction class.
Definition: GSkyDir.hpp:62
int dir2inx(const GSkyDir &dir) const
Returns pixel index for a given sky direction.
Definition: GSkyMap.cpp:1445
void close(void)
Close FITS file.
Definition: GFits.cpp:1342
bool contains(const GSkyDir &dir) const
Check if direction is contained in one of the regions.
#define G_CONSTRUCT_MAP
Definition: GSkyMap.cpp:56
GVector atan(const GVector &vector)
Computes arctan of vector elements.
Definition: GVector.cpp:1148
#define G_SHAPE
Definition: GSkyMap.cpp:59
Circular sky region class interface definition.
bool operator!=(const GEbounds &a, const GEbounds &b)
Energy boundaries inequality operator friend.
Definition: GEbounds.hpp:213
Abstract sky region base class interface definition.
GNdarray flux(void) const
Returns array with total flux for sky maps.
Definition: GSkyMap.cpp:1820
#define G_READ_WCS
Definition: GSkyMap.cpp:86
FITS table double column.
#define G_OP_ACCESS_2D
Definition: GSkyMap.cpp:66
Abstract world coordinate system base class definition.
virtual GSkyPixel dir2pix(const GSkyDir &dir) const =0
virtual GSkyProjection * clone(void) const =0
Clones object.
Mathematical function definitions.
#define G_READ_HEALPIX
Definition: GSkyMap.cpp:85
GVector log10(const GVector &vector)
Computes base10 logarithm of vector elements.
Definition: GVector.cpp:1295
Sky directions container class.
Definition: GSkyDirs.hpp:49
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489
FITS table abstract base class interface definition.