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