GammaLib  2.0.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-2020 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 Return sky region circle that encloses the sky map
2271  *
2272  * @return Enclosing sky region circle.
2273  *
2274  * Returns a sky region circle that encloses the sky map.
2275  *
2276  * For a sky map in HealPix projection the method returns a region circle
2277  * that encloses the full sky, with a centre at Right Ascension and
2278  * Declination of 0 and a radius of 180 degrees.
2279  *
2280  * For a sky map in World Coordinate projection the method returns a region
2281  * circle that is centred on the central pixel of the sky map (i.e. the
2282  * pixel with indices \f$(n_x/2,n_y/2)\f$, where \f$n_x\f$ and \f$n_y\f$ are
2283  * the number of pixels in the X- and Y-drection) and a radius that is
2284  * determined from the maximum radial distance of all pixels with respect
2285  * to the central pixel.
2286  ***************************************************************************/
2288 {
2289  // Initialise sky region circle
2290  GSkyRegionCircle region;
2291 
2292  // If the sky map is in HealPix projection then set radius to 180 deg
2293  if (projection()->code() == "HPX") {
2294  region = GSkyRegionCircle(0.0, 0.0, 180.0);
2295  }
2296 
2297  // ... otherwise sky map is in World Coordinate Projection
2298  else {
2299 
2300  // Initialise maximum radius
2301  double max_radius = 0.0;
2302  int max_index = -1;
2303 
2304  // Determine map centre
2305  GSkyPixel pixel(nx()/2.0, ny()/2.0);
2306  GSkyDir centre = pix2dir(pixel);
2307 
2308  // Determine maximum pixel distance with respect to map centre
2309  for (int i = 0; i < npix(); ++i) {
2310  double radius = inx2dir(i).dist_deg(centre);
2311  if (radius > max_radius) {
2312  max_radius = radius;
2313  max_index = i;
2314  }
2315  }
2316 
2317  // Consider edges of pixel with maximum distance
2318  if (max_index != -1) {
2319 
2320  // Set upper pixel edge just beyond the upper boundary since the
2321  // upper boundary is excluded
2322  const double up = 0.4999999;
2323 
2324  // Get pixel with maximum distance
2325  GSkyPixel pixel = inx2pix(max_index);
2326 
2327  // Get boundaries of pixel
2328  GSkyDir boundary1 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()-0.5));
2329  GSkyDir boundary2 = pix2dir(GSkyPixel(pixel.x(), pixel.y()-0.5));
2330  GSkyDir boundary3 = pix2dir(GSkyPixel(pixel.x()+up, pixel.y()-0.5));
2331  GSkyDir boundary4 = pix2dir(GSkyPixel(pixel.x()+up, pixel.y()));
2332  GSkyDir boundary5 = pix2dir(GSkyPixel(pixel.x()+up, pixel.y()+up));
2333  GSkyDir boundary6 = pix2dir(GSkyPixel(pixel.x(), pixel.y()+up));
2334  GSkyDir boundary7 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()+up));
2335  GSkyDir boundary8 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()));
2336 
2337  // Get radii
2338  double radius1 = boundary1.dist_deg(centre);
2339  double radius2 = boundary2.dist_deg(centre);
2340  double radius3 = boundary3.dist_deg(centre);
2341  double radius4 = boundary4.dist_deg(centre);
2342  double radius5 = boundary5.dist_deg(centre);
2343  double radius6 = boundary6.dist_deg(centre);
2344  double radius7 = boundary7.dist_deg(centre);
2345  double radius8 = boundary8.dist_deg(centre);
2346 
2347  // Determine maximum radius
2348  if (radius1 > max_radius) {
2349  max_radius = radius1;
2350  }
2351  if (radius2 > max_radius) {
2352  max_radius = radius2;
2353  }
2354  if (radius3 > max_radius) {
2355  max_radius = radius3;
2356  }
2357  if (radius4 > max_radius) {
2358  max_radius = radius4;
2359  }
2360  if (radius5 > max_radius) {
2361  max_radius = radius5;
2362  }
2363  if (radius6 > max_radius) {
2364  max_radius = radius6;
2365  }
2366  if (radius7 > max_radius) {
2367  max_radius = radius7;
2368  }
2369  if (radius8 > max_radius) {
2370  max_radius = radius8;
2371  }
2372 
2373  } // endif: considered edges of pixel
2374 
2375  // Set sky region
2376  region = GSkyRegionCircle(centre, max_radius);
2377 
2378  } // endelse: sky map was in World Coordinate Projection
2379 
2380  // Return region
2381  return region;
2382 }
2383 
2384 
2385 /***********************************************************************//**
2386  * @brief Load skymap from FITS file.
2387  *
2388  * @param[in] filename FITS file name..
2389  *
2390  * Loads HEALPix and non HEALPix skymaps. First searches for HEALPix map in
2391  * FITS file by scanning all HDUs for PIXTYPE=HEALPIX. If no HEALPix map has
2392  * been found then search load first non-empty image.
2393  *
2394  * @todo Do we have to restrict a HEALPix map to a BinTable and a WCS map
2395  * to a Double precision image???
2396  ***************************************************************************/
2397 void GSkyMap::load(const GFilename& filename)
2398 {
2399  // Clear sky map
2400  clear();
2401 
2402  // Initialize load flag
2403  bool loaded(false);
2404 
2405  // Open FITS file
2406  GFits fits(filename);
2407 
2408  // If an extension name is specified then first try loading that
2409  // extension
2410  if (filename.has_extname()) {
2411  const GFitsHDU& hdu = *fits.at(filename.extname());
2412  if (is_healpix(hdu)) {
2413  read_healpix(static_cast<const GFitsTable&>(hdu));
2414  loaded = true;
2415  }
2416  else if (is_wcs(hdu)) {
2417  read_wcs(static_cast<const GFitsImage&>(hdu));
2418  loaded = true;
2419  }
2420  }
2421 
2422  // ... otherwise is an extension number is specified then try loading
2423  // the corresponding HDU
2424  else if (filename.has_extno()) {
2425  const GFitsHDU& hdu = *fits.at(filename.extno());
2426  if (is_healpix(hdu)) {
2427  read_healpix(static_cast<const GFitsTable&>(hdu));
2428  loaded = true;
2429  }
2430  else if (is_wcs(hdu)) {
2431  read_wcs(static_cast<const GFitsImage&>(hdu));
2432  loaded = true;
2433  }
2434  }
2435 
2436  // If no map has yet been loaded then scan the file for an appropriate
2437  // HDU and read sky map from the first suitable HDU
2438  if (!loaded) {
2439 
2440  // Get number of HDUs
2441  int num = fits.size();
2442 
2443  // First search for HEALPix extension. We can skip the first
2444  // extension since this is always an image and a HEALPix map
2445  // is stored in a binary table
2446  for (int extno = 1; extno < num; ++extno) {
2447 
2448  // Get reference to HDU
2449  const GFitsHDU& hdu = *fits.at(extno);
2450 
2451  // If HDU is HEALPix then read data
2452  if (is_healpix(hdu)) {
2453  read_healpix(static_cast<const GFitsTable&>(hdu));
2454  loaded = true;
2455  break;
2456  }
2457 
2458  } // endfor: looped over HDUs
2459 
2460  // If we have not found a HEALPIX map then search now for an
2461  // image.
2462  if (!loaded) {
2463  for (int extno = 0; extno < num; ++extno) {
2464 
2465  // Get referene to HDU
2466  const GFitsHDU& hdu = *fits.at(extno);
2467 
2468  // If HDU is WCS then read data
2469  if (is_wcs(hdu)) {
2470  read_wcs(static_cast<const GFitsImage&>(hdu));
2471  //loaded = true;
2472  break;
2473  }
2474 
2475  } // endfor: looped over HDUs
2476  } // endif: no sky map yet loaded
2477  } // endif: no sky map yet loaded
2478 
2479  // Close FITS file
2480  fits.close();
2481 
2482  // Return
2483  return;
2484 }
2485 
2486 
2487 /***********************************************************************//**
2488  * @brief Save sky map into FITS file
2489  *
2490  * @param[in] filename FITS file name.
2491  * @param[in] clobber Overwrite existing file?
2492  *
2493  * Saves the sky map into a FITS file. If the file exists already and the
2494  * @p clobber parameter is true, the method will overwrite the content of
2495  * the existing file. Otherwise, an exception is thrown.
2496  ***************************************************************************/
2497 void GSkyMap::save(const GFilename& filename, const bool& clobber) const
2498 {
2499  // Continue only if we have data to save
2500  if (m_proj != NULL) {
2501 
2502  // Initialise HDU pointer
2503  GFitsHDU* hdu = NULL;
2504 
2505  // Case A: Skymap is Healpix
2506  if (m_proj->code() == "HPX") {
2507  hdu = create_healpix_hdu();
2508  }
2509 
2510  // Case B: Skymap is not Healpix
2511  else {
2512  hdu = create_wcs_hdu();
2513  }
2514 
2515  // If we have a valid HDU then save it now to the FITS file
2516  if (hdu != NULL) {
2517 
2518  // Set extension name
2519  if (filename.has_extname()) {
2520  hdu->extname(filename.extname());
2521  }
2522 
2523  // Initialise empty FITS file
2524  GFits fits;
2525 
2526  // Append sky map to FITS file
2527  fits.append(*hdu);
2528 
2529  // Save FITS file (without extension name which was extracted
2530  // earlier and set in the HDU)
2531  fits.saveto(filename.url(), clobber);
2532 
2533  // Delete HDU
2534  delete hdu;
2535 
2536  } // endif: HDU was valid
2537 
2538  } // endif: we had data to save
2539 
2540  // Return
2541  return;
2542 }
2543 
2544 
2545 /***********************************************************************//**
2546  * @brief Read skymap from FITS HDU
2547  *
2548  * @param[in] hdu FITS HDU.
2549  ***************************************************************************/
2550 void GSkyMap::read(const GFitsHDU& hdu)
2551 {
2552  // Free memory and initialise members
2553  free_members();
2554  init_members();
2555 
2556  // If HDU is a Healpix HDU then load a Healpix map
2557  if (is_healpix(hdu)) {
2558  read_healpix(static_cast<const GFitsTable&>(hdu));
2559  }
2560 
2561  // ... otherwise try loading as WCS map
2562  else if (is_wcs(hdu)) {
2563  read_wcs(static_cast<const GFitsImage&>(hdu));
2564  }
2565 
2566  // Return
2567  return;
2568 }
2569 
2570 
2571 /***********************************************************************//**
2572  * @brief Write sky map into FITS file
2573  *
2574  * @param[in] file FITS file pointer.
2575  * @param[in] extname Sky map extension name.
2576  * @return Pointer to written HDU.
2577  *
2578  * Write the sky map into a FITS file. Optionally, the extension name of
2579  * the FITS HDU can be specified using the @p extname parameter. The method
2580  * returns a pointer to the appended HDU. If no HDU has been appended, for
2581  * example because the sky map is empty, the method will return NULL.
2582  ***************************************************************************/
2583 GFitsHDU* GSkyMap::write(GFits& file, const std::string& extname) const
2584 {
2585  // Initialise pointer to appended HDU
2586  GFitsHDU* hdu_appended = NULL;
2587 
2588  // Continue only if we have data to save
2589  if (m_proj != NULL) {
2590 
2591  // Initialise local HDU pointer
2592  GFitsHDU* hdu = NULL;
2593 
2594  // Case A: Skymap is Healpix
2595  if (m_proj->code() == "HPX") {
2596  hdu = create_healpix_hdu();
2597  }
2598 
2599  // Case B: Skymap is not Healpix
2600  else {
2601  hdu = create_wcs_hdu();
2602  }
2603 
2604  // Append HDU to FITS file.
2605  if (hdu != NULL) {
2606 
2607  // Optionally set extension name
2608  if (extname.length() > 0) {
2609  hdu->extname(extname);
2610  }
2611 
2612  // Append HDU and recover the pointer to the appended HDU
2613  hdu_appended = file.append(*hdu);
2614 
2615  // Delete HDU
2616  delete hdu;
2617 
2618  } // endif: there was a valid HDU
2619 
2620  } // endif: we had data to save
2621 
2622  // Return pointer to appended HDU
2623  return hdu_appended;
2624 }
2625 
2626 
2627 /***********************************************************************//**
2628  * @brief Publish sky map
2629  *
2630  * @param[in] name Name of sky map.
2631  *
2632  * Publishes the sky map on a Virtual Observatory Hub. If no Hub is currently
2633  * active, the method will start a new Hub. If on sky map has been allocated
2634  * the method does nothing.
2635  ***************************************************************************/
2636 void GSkyMap::publish(const std::string& name) const
2637 {
2638  // Create FITS file containing the sky map
2639  GFits fits;
2640 
2641  // Write sky map into FITS file and return a pointer to the written HDU.
2642  // This method returns a NULL pointer if no sky map projection exists,
2643  // which typically is true for empty sky maps.
2644  GFitsHDU* hdu = write(fits);
2645 
2646  // Continue only if HDU pointer is valid
2647  if (hdu != NULL) {
2648 
2649  // Optionally set extension name
2650  if (!name.empty()) {
2651  hdu->extname(name);
2652  }
2653 
2654  // Create VO Client
2655  GVOClient client;
2656 
2657  // Publish map using VO client
2658  client.publish(*hdu);
2659 
2660  } // endif: HDU pointer was valid
2661 
2662  // Return
2663  return;
2664 }
2665 
2666 
2667 /***********************************************************************//**
2668  * @brief Print sky map
2669  *
2670  * @param[in] chatter Chattiness.
2671  * @return String containing sky map information.
2672  ***************************************************************************/
2673 std::string GSkyMap::print(const GChatter& chatter) const
2674 {
2675  // Initialise result string
2676  std::string result;
2677 
2678  // Continue only if chatter is not silent
2679  if (chatter != SILENT) {
2680 
2681  // Append header
2682  result.append("=== GSkyMap ===");
2683 
2684  // Append number of pixels and number of maps
2685  result.append("\n"+gammalib::parformat("Number of pixels"));
2686  result.append(gammalib::str(m_num_pixels));
2687 
2688  // Append WCS dimension information
2689  if (m_proj != NULL && m_proj->size() == 2) {
2690  result.append("\n"+gammalib::parformat("X axis dimension"));
2691  result.append(gammalib::str(m_num_x));
2692  result.append("\n"+gammalib::parformat("Y axis dimension"));
2693  result.append(gammalib::str(m_num_y));
2694  }
2695 
2696  // Append number of maps
2697  result.append("\n"+gammalib::parformat("Number of maps"));
2698  result.append(gammalib::str(m_num_maps));
2699 
2700  // Append shape of maps
2701  result.append("\n"+gammalib::parformat("Shape of maps"));
2702  if (m_shape.size() > 0) {
2703  std::string shape = "(";
2704  for (int i = 0; i < m_shape.size(); ++i) {
2705  if (i > 0) {
2706  shape += ", ";
2707  }
2708  shape += gammalib::str(m_shape[i]);
2709  }
2710  shape += ")";
2711  result.append(shape);
2712  }
2713  else {
2714  result.append("not defined");
2715  }
2716 
2717  // Append sky projection information
2718  if (gammalib::reduce(chatter) > SILENT) {
2719  if (m_proj != NULL) {
2720  result.append("\n"+m_proj->print(gammalib::reduce(chatter)));
2721  }
2722  else {
2723  result.append("\n"+gammalib::parformat("Sky projection"));
2724  result.append("not defined");
2725  }
2726  }
2727 
2728  } // endif: chatter was not silent
2729 
2730  // Return result
2731  return result;
2732 }
2733 
2734 /*==========================================================================
2735  = =
2736  = Private methods =
2737  = =
2738  ==========================================================================*/
2739 
2740 /***********************************************************************//**
2741  * @brief Initialise class members
2742  ***************************************************************************/
2744 {
2745  // Initialise members
2746  m_num_pixels = 0;
2747  m_num_maps = 0;
2748  m_num_x = 0;
2749  m_num_y = 0;
2750  m_shape.clear();
2751  m_proj = NULL;
2752  m_pixels.clear();
2753 
2754  // Initialise computation cache
2755  m_hascache = false;
2756  m_contained = false;
2757  m_last_dir.clear();
2758  m_interpol.clear();
2759 
2760  // Return
2761  return;
2762 }
2763 
2764 
2765 /***********************************************************************//**
2766  * @brief Copy class members
2767  *
2768  * @param[in] map Sky map.
2769  ***************************************************************************/
2771 {
2772  // Copy attributes
2773  m_num_pixels = map.m_num_pixels;
2774  m_num_maps = map.m_num_maps;
2775  m_num_x = map.m_num_x;
2776  m_num_y = map.m_num_y;
2777  m_shape = map.m_shape;
2778  m_pixels = map.m_pixels;
2779 
2780  // Copy computation cache
2781  m_hascache = map.m_hascache;
2782  m_contained = map.m_contained;
2783  m_last_dir = map.m_last_dir;
2784  m_interpol = map.m_interpol;
2785 
2786  // Clone sky projection if it is valid
2787  if (map.m_proj != NULL) m_proj = map.m_proj->clone();
2788 
2789  // Return
2790  return;
2791 }
2792 
2793 
2794 /***********************************************************************//**
2795  * @brief Delete class members
2796  ***************************************************************************/
2798 {
2799  // Free memory
2800  if (m_proj != NULL) {
2801  delete m_proj;
2802  }
2803 
2804  // Signal free pointers
2805  m_proj = NULL;
2806 
2807  // Return
2808  return;
2809 }
2810 
2811 
2812 /***********************************************************************//**
2813  * @brief Set World Coordinate System
2814  *
2815  * @param[in] wcs World Coordinate System code.
2816  * @param[in] coords Coordinate system.
2817  * @param[in] crval1 X value of reference pixel.
2818  * @param[in] crval2 Y value of reference pixel.
2819  * @param[in] crpix1 X index of reference pixel.
2820  * @param[in] crpix2 Y index of reference pixel.
2821  * @param[in] cdelt1 Increment in x direction at reference pixel (deg).
2822  * @param[in] cdelt2 Increment in y direction at reference pixel (deg).
2823  * @param[in] cd Astrometry parameters (2x2 matrix, deg/pixel).
2824  * @param[in] pv2 Projection parameters (length WCS type dependent).
2825  *
2826  * @exception GException::wcs_invalid
2827  * Invalid wcs parameter (World Coordinate System not known).
2828  *
2829  * This method sets the WCS projection pointer based on the WCS code and
2830  * sky map parameters. It makes use of the GWcsRegistry class to allocate
2831  * the correct derived class. Note that this method does not support the
2832  * HPX projection.
2833  *
2834  * @todo Remove cd and pv2 parameters.
2835  ***************************************************************************/
2836 void GSkyMap::set_wcs(const std::string& wcs, const std::string& coords,
2837  const double& crval1, const double& crval2,
2838  const double& crpix1, const double& crpix2,
2839  const double& cdelt1, const double& cdelt2,
2840  const GMatrix& cd, const GVector& pv2)
2841 {
2842  // Convert WCS to upper case
2843  std::string uwcs = gammalib::toupper(wcs);
2844 
2845  // Check if HPX was requested (since this is not allowed)
2846  if (uwcs == "HPX") {
2847  throw GException::wcs_invalid(G_SET_WCS, uwcs,
2848  "Method not valid for HPX projection.");
2849  }
2850 
2851  // ... otherwise get projection from registry
2852  else {
2853  // Allocate WCS registry
2854  GWcsRegistry registry;
2855 
2856  // Allocate projection from registry
2857  GWcs* projection = registry.alloc(uwcs);
2858  m_proj = projection;
2859 
2860  // Signal if projection type is not known
2861  if (projection == NULL) {
2862  std::string message = "Projection code not known. "
2863  "Should be one of "+registry.list()+".";
2864  throw GException::wcs_invalid(G_SET_WCS, uwcs, message);
2865  }
2866 
2867  // Setup WCS
2868  projection->set(coords, crval1, crval2, crpix1, crpix2,
2869  cdelt1, cdelt2);
2870 
2871  } // endelse: got projection from registry
2872 
2873  // Return
2874  return;
2875 }
2876 
2877 
2878 /***********************************************************************//**
2879  * @brief Read Healpix data from FITS table.
2880  *
2881  * @param[in] table FITS table.
2882  *
2883  * HEALPix data may be stored in various formats depending on the
2884  * application that has writted the data. HEALPix IDL, for example, may
2885  * store the data in vectors of length 1024 if the number of pixels is
2886  * a multiple of 1024. On the other hand, vectors may also be used to store
2887  * several HEALPix maps into a single column. Alternatively, multiple maps
2888  * may be stored in multiple columns.
2889  ***************************************************************************/
2891 {
2892  // Determine number of rows and columns in table
2893  int nrows = table.nrows();
2894  int ncols = table.ncols();
2895  #if defined(G_READ_HEALPIX_DEBUG)
2896  std::cout << "nrows=" << nrows << " ncols=" << ncols << std::endl;
2897  #endif
2898 
2899  // Allocate Healpix projection
2900  m_proj = new GHealpix;
2901 
2902  // Read WCS information from FITS header
2903  m_proj->read(table);
2904 
2905  // Set number of pixels based on NSIDE parameter
2906  m_num_pixels = static_cast<GHealpix*>(m_proj)->npix();
2907  #if defined(G_READ_HEALPIX_DEBUG)
2908  std::cout << "m_num_pixels=" << m_num_pixels << std::endl;
2909  #endif
2910 
2911  // Number of map pixels has to be a multiple of the number of
2912  // rows in column
2913  if (m_num_pixels % nrows != 0) {
2915  m_num_pixels);
2916  }
2917 
2918  // Determine vector length for HEALPix data storage
2919  int nentry = m_num_pixels / nrows;
2920  #if defined(G_READ_HEALPIX_DEBUG)
2921  std::cout << "nentry=" << nentry << std::endl;
2922  #endif
2923 
2924  // Determine number of maps from the number of maps that fit into
2925  // all columns. Only count columns that can fully hold the map.
2926  m_num_maps = 0;
2927  for (int icol = 0; icol < ncols; ++icol) {
2928  const GFitsTableCol* col = table[icol];
2929  if (col->number() % nentry == 0) {
2930  m_num_maps += col->number() / nentry;
2931  }
2932  }
2933  #if defined(G_READ_HEALPIX_DEBUG)
2934  std::cout << "m_num_maps=" << m_num_maps << std::endl;
2935  #endif
2936 
2937  // Initialise shape of maps
2938  m_shape.clear();
2939  m_shape.push_back(m_num_maps);
2940 
2941  // Allocate pixels to hold the map
2943 
2944  // Initialise map counter
2945  int imap = 0;
2946 
2947  // Loop over all columns
2948  for (int icol = 0; icol < ncols; ++icol) {
2949 
2950  // Get next column
2951  const GFitsTableCol* col = table[icol];
2952 
2953  // Only consider columns that can fully hold maps
2954  if (col->number() % nentry == 0) {
2955 
2956  // Determine number of maps in column
2957  int num = col->number() / nentry;
2958 
2959  // Loop over all maps in column
2960  int inx_start = 0;
2961  int inx_end = nentry;
2962  for (int i = 0; i < num; ++i) {
2963 
2964  // Load map
2965  double *ptr = m_pixels.data() + m_num_pixels*imap;
2966  for (int row = 0; row < col->nrows(); ++row) {
2967  for (int inx = inx_start; inx < inx_end; ++inx) {
2968  *ptr++ = col->real(row,inx);
2969  }
2970  }
2971  #if defined(G_READ_HEALPIX_DEBUG)
2972  std::cout << "Load map=" << imap << " index="
2973  << inx_start << "-" << inx_end << std::endl;
2974  #endif
2975 
2976  // Increment index range
2977  inx_start = inx_end;
2978  inx_end += nentry;
2979 
2980  // Increment map counter
2981  imap++;
2982 
2983  // Break if we have loaded all maps
2984  if (imap >= m_num_maps) {
2985  break;
2986  }
2987 
2988  } // endfor: looped over all maps in column
2989  } // endif: column could fully hold maps
2990 
2991  // Break if we have loaded all maps
2992  if (imap >= m_num_maps) {
2993  break;
2994  }
2995 
2996  } // endfor: looped over all columns
2997 
2998  // Return
2999  return;
3000 }
3001 
3002 
3003 /***********************************************************************//**
3004  * @brief Read WCS image from FITS HDU
3005  *
3006  * @param[in] image FITS image.
3007  *
3008  * @exception GException::invalid_argument
3009  * FITS image has less than two dimensions
3010  * @exception GException::invalid_value
3011  * Sky map covers more than 360 deg in longitude or 180 deg in
3012  * latitude
3013  *
3014  * Reads sky maps from a FITS image extension containing a set of maps given
3015  * in World Coordinate System. The method handles general n-dimensional
3016  * images and sets the map shape attribute according to the number of map
3017  * dimensions found in the FITS HDU.
3018  ***************************************************************************/
3019 void GSkyMap::read_wcs(const GFitsImage& image)
3020 {
3021  // Throw an exception if the FITS image is not at least a 2D image
3022  if (image.naxis() < 2) {
3023  std::string msg = "Sky map has "+gammalib::str(image.naxis())+
3024  " dimensions, which is less than the two dimensions"
3025  " that are required for a WCS image.";
3027  }
3028 
3029  // Allocate WCS
3030  alloc_wcs(image);
3031 
3032  // Read projection information from FITS header
3033  m_proj->read(image);
3034 
3035  // Clear map shape
3036  m_shape.clear();
3037 
3038  // Extract map dimension, number of maps and map shape from FITS image
3039  m_num_x = image.naxes(0);
3040  m_num_y = image.naxes(1);
3041  if (image.naxis() == 2) {
3042  m_num_maps = 1;
3043  m_shape.push_back(m_num_maps);
3044  }
3045  else {
3046  m_num_maps = image.naxes(2);
3047  m_shape.push_back(m_num_maps);
3048  for (int i = 3; i < image.naxis(); ++i) {
3049  m_num_maps *= image.naxes(i);
3050  m_shape.push_back(image.naxes(i));
3051  }
3052  }
3053  #if defined(G_READ_WCS_DEBUG)
3054  std::cout << "m_num_x=" << m_num_x << std::endl;
3055  std::cout << "m_num_y=" << m_num_y << std::endl;
3056  std::cout << "m_num_maps=" << m_num_maps << std::endl;
3057  #endif
3058 
3059  // Compute number of pixels
3061  #if defined(G_READ_WCS_DEBUG)
3062  std::cout << "m_num_pixels=" << m_num_pixels << std::endl;
3063  #endif
3064 
3065  // Allocate pixels to hold the map
3066  m_pixels = GNdarray(m_num_x, m_num_y, m_num_maps);
3067 
3068  // Read image
3069  double* ptr = m_pixels.data();
3070  int size = m_num_pixels * m_num_maps;
3071  for (int i = 0; i < size; ++i) {
3072  *ptr++ = (double)image.pixel(i);
3073  }
3074 
3075  // Check if the map is too large
3076  double x_size = std::abs(m_num_x * static_cast<GWcs*>(m_proj)->cdelt(0));
3077  double y_size = std::abs(m_num_y * static_cast<GWcs*>(m_proj)->cdelt(1));
3078  if (x_size > 360.001) {
3079  std::string msg = "Skymap covers "+gammalib::str(x_size)+" degrees "
3080  "in the X axis which results in ambiguous "
3081  "coordinate transformations. Please provide a "
3082  "skymap that does not cover more than 360 degrees.";
3084  }
3085  if (y_size > 180.001) {
3086  std::string msg = "Skymap covers "+gammalib::str(y_size)+" degrees "
3087  "in the Y axis which results in ambiguous "
3088  "coordinate transformations. Please provide a "
3089  "skymap that does not cover more than 180 degrees.";
3091  }
3092 
3093  // Return
3094  return;
3095 }
3096 
3097 
3098 /***********************************************************************//**
3099  * @brief Allocate WCS class
3100  *
3101  * @param[in] image FITS image.
3102  *
3103  * @exception GException::fits_key_not_found
3104  * Unable to find required FITS header keyword.
3105  * @exception GException::skymap_bad_ctype
3106  * CTYPE1 and CTYPE2 keywords are incompatible.
3107  * @exception GException::wcs_invalid
3108  * WCS projection of FITS file not supported by GammaLib.
3109  ***************************************************************************/
3110 void GSkyMap::alloc_wcs(const GFitsImage& image)
3111 {
3112  // Get standard keywords
3113  std::string ctype1 = image.string("CTYPE1");
3114  std::string ctype2 = image.string("CTYPE2");
3115 
3116  // Extract projection type
3117  std::string xproj = ctype1.substr(5,3);
3118  std::string yproj = ctype2.substr(5,3);
3119 
3120  // Check that projection type is identical on both axes
3121  if (xproj != yproj) {
3123  ctype1, ctype2);
3124  }
3125 
3126  // Allocate WCS registry
3127  GWcsRegistry registry;
3128 
3129  // Allocate projection from registry
3130  m_proj = registry.alloc(xproj);
3131 
3132  // Signal if projection type is not known
3133  if (m_proj == NULL) {
3134  std::string message = "Projection code not known. "
3135  "Should be one of "+registry.list()+".";
3136  throw GException::wcs_invalid(G_ALLOC_WCS, xproj, message);
3137  }
3138 
3139  // Return
3140  return;
3141 }
3142 
3143 
3144 /***********************************************************************//**
3145  * @brief Create FITS HDU containing Healpix data
3146  *
3147  * This method allocates a binary table HDU that contains the Healpix data.
3148  * Deallocation of the table has to be done by the client.
3149  ***************************************************************************/
3151 {
3152  // Initialise result to NULL pointer
3153  GFitsBinTable* hdu = NULL;
3154 
3155  // Compute size of Healpix data
3156  int size = m_num_pixels * m_num_maps;
3157 
3158  // Continue only if we have pixels
3159  if (size > 0) {
3160 
3161  // Set number of rows and columns
3162  int rows = m_num_pixels;
3163  int number = m_num_maps;
3164 
3165  // Create column to hold Healpix data
3166  GFitsTableDoubleCol column = GFitsTableDoubleCol("DATA", rows, number);
3167 
3168  // Fill data into column
3169  const double* ptr = m_pixels.data();
3170  for (int inx = 0; inx < number; ++inx) {
3171  for (int row = 0; row < rows; ++row) {
3172  column(row,inx) = *ptr++;
3173  }
3174  }
3175 
3176  // Create HDU that contains Healpix map in a binary table
3177  hdu = new GFitsBinTable(rows);
3178  hdu->append(column);
3179 
3180  } // endif: there were pixels
3181 
3182  // ... otherwise create an empty header
3183  else {
3184  hdu = new GFitsBinTable;
3185  }
3186 
3187  // Set extension name
3188  hdu->extname("HEALPIX");
3189 
3190  // If we have WCS information then write into FITS header
3191  if (m_proj != NULL) {
3192  m_proj->write(*hdu);
3193  }
3194 
3195  // Set additional keywords
3196  hdu->card("NBRBINS", m_num_maps, "Number of HEALPix maps");
3197 
3198  // Return HDU
3199  return hdu;
3200 }
3201 
3202 
3203 /***********************************************************************//**
3204  * @brief Create FITS HDU containing WCS image
3205  *
3206  * This method allocates an image HDU that contains the WCS image data.
3207  * Deallocation of the image has to be done by the client.
3208  *
3209  * @todo Set additional keywords.
3210  ***************************************************************************/
3212 {
3213  // Initialise result to NULL pointer
3214  GFitsImageDouble* hdu = NULL;
3215 
3216  // Compute size of Healpix data
3217  int size = m_num_pixels * m_num_maps;
3218 
3219  // Continue only if we have pixels
3220  if (size > 0) {
3221 
3222  // Set dimension vector of all axes. In case that only one map
3223  // exists then create simply a 2D image
3224  std::vector<int> naxes;
3225  naxes.push_back(m_num_x);
3226  naxes.push_back(m_num_y);
3227  if (m_num_maps > 1) {
3228  for (int i = 0; i < ndim(); ++i) {
3229  naxes.push_back(m_shape[i]);
3230  }
3231  }
3232 
3233  // Allocate image
3234  hdu = new GFitsImageDouble(naxes);
3235 
3236  // Store data in image
3237  if (naxes.size() == 2) {
3238  const double* ptr = m_pixels.data();
3239  for (int iy = 0; iy < m_num_y; ++iy) {
3240  for (int ix = 0; ix < m_num_x; ++ix) {
3241  (*hdu)(ix,iy) = *ptr++;
3242  }
3243  }
3244  }
3245  else {
3246  const double* ptr = m_pixels.data();
3247  for (int imap = 0; imap < m_num_maps; ++imap) {
3248  for (int iy = 0; iy < m_num_y; ++iy) {
3249  for (int ix = 0; ix < m_num_x; ++ix) {
3250  (*hdu)(ix,iy,imap) = *ptr++;
3251  }
3252  }
3253  }
3254  }
3255 
3256  } // endif: there were pixels
3257 
3258  // ... otherwise create an empty header
3259  else {
3260  hdu = new GFitsImageDouble;
3261  }
3262 
3263  // Set extension name
3264  hdu->extname("IMAGE");
3265 
3266  // If we have sky projection information then write into FITS header
3267  if (m_proj != NULL) {
3268  m_proj->write(*hdu);
3269  }
3270 
3271  // Set additional keywords
3272  //TODO
3273 
3274  // Return HDU
3275  return hdu;
3276 }
3277 
3278 
3279 /***********************************************************************//**
3280  * @brief Compute solid angle subtended by 4 sky directions
3281  *
3282  * @param[in] dir1 First sky direction.
3283  * @param[in] dir2 Second sky direction.
3284  * @param[in] dir3 Third sky direction.
3285  * @param[in] dir4 Forth sky direction.
3286  * @return Solid angle (steradians).
3287  *
3288  * Estimate the solid angle subtended by 4 sky directions using Huilier's
3289  * theorem.
3290  *
3291  * Below, the definiton of the pixel cornes and sides are shown as used
3292  * within the code.
3293  *
3294  * a12
3295  * 1---------2
3296  * |\ /|
3297  * | \a13 / |
3298  * | \ / |
3299  * | \ / |
3300  * a14| X |a23
3301  * | / \ |
3302  * | / \ |
3303  * | /a24 \ |
3304  * |/ \|
3305  * 4---------3
3306  * a34
3307  *
3308  ***************************************************************************/
3309 double GSkyMap::solidangle(const GSkyDir& dir1, const GSkyDir& dir2,
3310  const GSkyDir& dir3, const GSkyDir& dir4) const
3311 {
3312  // Initialise solid angle
3313  double solidangle = 0.0;
3314 
3315  // Compute angular distances between pixel corners
3316  double a12 = dir1.dist(dir2);
3317  double a14 = dir1.dist(dir4);
3318  double a23 = dir2.dist(dir3);
3319  double a24 = dir2.dist(dir4);
3320  double a34 = dir3.dist(dir4);
3321 
3322  // Special case: a12 or a14 is zero, then pixel is a triangle composed
3323  // of [2,3,4]
3324  if (a12 <= 0.0 || a14 <= 0.0) {
3325  double s = 0.5 * (a23 + a34 + a24);
3326  solidangle = 4.0 * std::atan(std::sqrt(std::tan(0.5*s) *
3327  std::tan(0.5*(s-a23)) *
3328  std::tan(0.5*(s-a34)) *
3329  std::tan(0.5*(s-a24))));
3330  }
3331 
3332  // Special case: a23 or a 34 is zero, then pixel is a triangle composed
3333  // of [1,2,4]
3334  else if (a23 <= 0.0 || a34 <= 0.0) {
3335  double s = 0.5 * (a12 + a24 + a14);
3336  solidangle = 4.0 * std::atan(std::sqrt(std::tan(0.5*s) *
3337  std::tan(0.5*(s-a12)) *
3338  std::tan(0.5*(s-a24)) *
3339  std::tan(0.5*(s-a14))));
3340  }
3341 
3342  // Otherwise we have a polygon
3343  else {
3344 
3345  // Triangle 1 [1,2,4]
3346  double s1 = 0.5 * (a12 + a24 + a14);
3347  double excess1 = std::atan(std::sqrt(std::tan(0.5*s1) *
3348  std::tan(0.5*(s1-a12)) *
3349  std::tan(0.5*(s1-a24)) *
3350  std::tan(0.5*(s1-a14))));
3351 
3352  // Triangle 2 [2,3,4]
3353  double s2 = 0.5 * (a23 + a34 + a24);
3354  double excess2 = std::atan(std::sqrt(std::tan(0.5*s2) *
3355  std::tan(0.5*(s2-a23)) *
3356  std::tan(0.5*(s2-a34)) *
3357  std::tan(0.5*(s2-a24))));
3358 
3359  // Determine solid angle
3360  solidangle = 4.0 * (excess1 + excess2);
3361 
3362  } // endif: we had a polynom
3363 
3364  // Return solid angle
3365  return solidangle;
3366 }
3367 
3368 
3369 /***********************************************************************//**
3370  * @brief Compute solid angle subtended by 3 sky directions
3371  *
3372  * @param[in] dir1 First sky direction.
3373  * @param[in] dir2 Second sky direction.
3374  * @param[in] dir3 Third sky direction.
3375  * @return Solid angle (steradians).
3376  *
3377  * Estimate the solid angle subtended by 3 sky directions using Huilier's
3378  * theorem.
3379  *
3380  * Below, the definiton of the pixel cornes and sides are shown as used
3381  * within the code.
3382  *
3383  * a12
3384  * 1---------2
3385  * | /
3386  * | /
3387  * | /
3388  * | /
3389  * a13| /a23
3390  * | /
3391  * | /
3392  * | /
3393  * |/
3394  * 3
3395  *
3396  ***************************************************************************/
3397 double GSkyMap::solidangle(const GSkyDir& dir1, const GSkyDir& dir2,
3398  const GSkyDir& dir3) const
3399 {
3400  // Compute angular distances between pixel corners
3401  double a12 = dir1.dist(dir2);
3402  double a13 = dir1.dist(dir3);
3403  double a23 = dir2.dist(dir3);
3404 
3405  // Compute solid angle
3406  double s = 0.5 * (a12 + a23 + a13);
3407  double solidangle = 4.0 * std::atan(std::sqrt(std::tan(0.5*s) *
3408  std::tan(0.5*(s-a12)) *
3409  std::tan(0.5*(s-a23)) *
3410  std::tan(0.5*(s-a13))));
3411 
3412  // Return solid angle
3413  return solidangle;
3414 }
3415 
3416 
3417 /***********************************************************************//**
3418  * @brief Checks whether a circular region overlaps with this map
3419  *
3420  * @param[in] region Circular region
3421  * @return True if circular region overlaps with map, false otherwise
3422  *
3423  * The check is done by first testing whether the central pointing
3424  * position of the observation falls into the sky map. If this is false,
3425  * then pixel positions that border the sky map are tested for whether or
3426  * not they fall into the observation region. The positions tested can be
3427  * visualized as follows, where '*' marks the positions tested.
3428  *
3429  * * * * * * *
3430  * +---+---+---+---+
3431  * * |0,2|1,2|2,2|3,2| *
3432  * +---+---+---+---+
3433  * * |0,1|1,1|2,1|3,1| *
3434  * +---+---+---+---+
3435  * * |0,0|1,0|2,0|3,0| *
3436  * +---+---+---+---+
3437  * * * * * * *
3438  *
3439  ***************************************************************************/
3441 {
3442  // Initialise overlap with true
3443  bool overlap = false;
3444 
3445  // If the centre of the region is inside the map then signal overlap ...
3446  if (contains(region.centre())) {
3447  overlap = true;
3448  }
3449 
3450  // ... otherwise loop through each of the outer bins and check if the
3451  // region overlaps with them (within some delta)
3452  else {
3453 
3454  // Initialise test pixels
3455  GSkyPixel pix1(0.0,0.0);
3456  GSkyPixel pix2(0.0,0.0);
3457 
3458  // Check the bottom and top bins
3459  pix1.y(-1.0);
3460  pix2.y(double(ny()));
3461  for (double xbin = -1.0; xbin <= nx(); xbin += 1.0) {
3462  pix1.x(xbin);
3463  pix2.x(xbin);
3464  if ((region.contains(pix2dir(pix1))) ||
3465  (region.contains(pix2dir(pix2)))) {
3466  overlap = true;
3467  break;
3468  }
3469  } // endfor: loop on xbin
3470 
3471  // If no overlap has been found then check now the left and right
3472  // bins
3473  if (!overlap) {
3474  pix1.x(-1.0);
3475  pix2.x(double(nx()));
3476  for (double ybin = 0.0; ybin < ny(); ybin += 1.0) {
3477  pix1.y(ybin);
3478  pix2.y(ybin);
3479  if ((region.contains(pix2dir(pix1))) ||
3480  (region.contains(pix2dir(pix2)))) {
3481  overlap = true;
3482  break;
3483  }
3484  } // endfor: loop on ybin
3485  }
3486 
3487  } // endelse: check border bins
3488 
3489  // Return overlap flag
3490  return overlap;
3491 }
3492 
3493 
3494 /***********************************************************************//**
3495  * @brief Check if HDU contains HEALPix data
3496  *
3497  * @param[in] hdu FITS Header Data Unit.
3498  * @return True is HDU contains HEALPix data.
3499  *
3500  * Returns true if the HDU is not an image and the HDU has the "PIXTYPE"
3501  * keyword set to "HEALPIX".
3502  ***************************************************************************/
3503 bool GSkyMap::is_healpix(const GFitsHDU& hdu) const
3504 {
3505  // Initialise flag
3506  bool flag(false);
3507 
3508  // If PIXTYPE keyword equals "HEALPIX" then signal that we have
3509  // HEALPix data
3510  if ((hdu.exttype() != GFitsHDU::HT_IMAGE) &&
3511  (hdu.has_card("PIXTYPE") && (hdu.string("PIXTYPE") == "HEALPIX"))) {
3512  flag = true;
3513  }
3514 
3515  // Return flag
3516  return (flag);
3517 }
3518 
3519 
3520 /***********************************************************************//**
3521  * @brief Check if HDU contains WCS data
3522  *
3523  * @param[in] hdu FITS Header Data Unit.
3524  * @return True is HDU contains WCS data.
3525  *
3526  * Returns true if the HDU is an image and the HDU has the "NAXIS" keyword
3527  * set to a value >= 2.
3528  ***************************************************************************/
3529 bool GSkyMap::is_wcs(const GFitsHDU& hdu) const
3530 {
3531  // Initialise flag
3532  bool flag(false);
3533 
3534  // If extension is an image thn signal that we have WCS data
3535  if ((hdu.exttype() == GFitsHDU::HT_IMAGE) &&
3536  (hdu.has_card("NAXIS") && (hdu.integer("NAXIS") >= 2))) {
3537  flag = true;
3538  }
3539 
3540  // Return flag
3541  return (flag);
3542 }
3543 
3544 
3545 /***********************************************************************//**
3546  * @brief Check if map is the same
3547  *
3548  * @param[in] map Sky map.
3549  * @return True is sky map definition is identical.
3550  *
3551  * Returns true if the sky map definition is identical to the one of the
3552  * current map.
3553  ***************************************************************************/
3554 bool GSkyMap::is_same(const GSkyMap& map) const
3555 {
3556  // Initialise flag
3557  bool flag = true;
3558 
3559  // Single loop that can be exited any time
3560  do {
3561 
3562  // If sky projection differs then set flag to false and break
3563  if (m_proj->code() != map.m_proj->code()) {
3564  flag = false;
3565  break;
3566  }
3567 
3568  // If sky projection has different coordinate system then set flag
3569  // to false and break
3570  if (m_proj->coordsys() != map.m_proj->coordsys()) {
3571  flag = false;
3572  break;
3573  }
3574 
3575  // If sky map dimension differs then set flag to false and break
3576  if ((m_num_x != map.m_num_x) || (m_num_y != map.m_num_y)) {
3577  flag = false;
3578  break;
3579  }
3580 
3581  // For WCS projections, if sky maps have different reference values,
3582  // pixels or pixel sizes then set flag to false and break
3583  const GWcs* wcs1 = dynamic_cast<const GWcs*>(m_proj);
3584  const GWcs* wcs2 = dynamic_cast<const GWcs*>(map.m_proj);
3585  if ((wcs1 != NULL) && (wcs2 != NULL)) {
3586  if ((wcs1->crval(0) != wcs2->crval(0)) ||
3587  (wcs1->crval(1) != wcs2->crval(1)) ||
3588  (wcs1->crpix(0) != wcs2->crpix(0)) ||
3589  (wcs1->crpix(1) != wcs2->crpix(1)) ||
3590  (wcs1->cdelt(0) != wcs2->cdelt(0)) ||
3591  (wcs1->cdelt(1) != wcs2->cdelt(1))) {
3592  flag = false;
3593  break;
3594  }
3595  }
3596 
3597  // For Healpix projections, if sky maps have different number of pixels,
3598  // nside or ordering then set flag to false and break
3599  const GHealpix* healpix1 = dynamic_cast<const GHealpix*>(m_proj);
3600  const GHealpix* healpix2 = dynamic_cast<const GHealpix*>(map.m_proj);
3601  if ((healpix1 != NULL) && (healpix2 != NULL)) {
3602  if ((healpix1->npix() != healpix2->npix()) ||
3603  (healpix1->nside() != healpix2->nside()) ||
3604  (healpix1->ordering() != healpix2->ordering())) {
3605  flag = false;
3606  break;
3607  }
3608  }
3609 
3610  } while(false);
3611 
3612  // Return flag
3613  return (flag);
3614 }
3615 
3616 
3617 /***********************************************************************//**
3618  * @brief Return smoothing kernel
3619  *
3620  * @param[in] kernel Smoothing kernel type ("DISK", "GAUSSIAN").
3621  * @param[in] par Smoothing parameter (>0).
3622  * @return Array filled with smoothing kernel.
3623  ***************************************************************************/
3624 GNdarray GSkyMap::smooth_kernel(const std::string& kernel,
3625  const double& par) const
3626 {
3627  // Get pointer on WCS projection
3628  const GWcs* wcs = dynamic_cast<const GWcs*>(m_proj);
3629  if (wcs == NULL) {
3630  std::string msg = "Sky map is not a WCS projection. Method is only "
3631  "valid for WCS projections.";
3633  }
3634 
3635  // Get X and Y step size
3636  double dx = wcs->cdelt(0);
3637  double dy = wcs->cdelt(1);
3638 
3639  // Initialise kernel
3640  GNdarray kern(m_num_x, m_num_y);
3641 
3642  // Initialise kernel sum
3643  double sum = 0.0;
3644 
3645  // Handle disk kernel
3646  if (gammalib::toupper(kernel) == "DISK") {
3647 
3648  // Fill kernel
3649  for (int ix1 = 0, ix2 = m_num_x; ix1 < m_num_x; ++ix1, --ix2) {
3650  double x = ix1 * dx;
3651  double xqs = x * x;
3652  for (int iy1 = 0, iy2 = m_num_y; iy1 < m_num_y; ++iy1, --iy2) {
3653  double y = iy1 * dy;
3654  if (std::sqrt(xqs + y*y) <= par) {
3655  kern(ix1,iy1) += 1.0;
3656  sum += 1.0;
3657  if (ix2 < m_num_x) {
3658  kern(ix2,iy1) += 1.0;
3659  sum += 1.0;
3660  }
3661  if (iy2 < m_num_y) {
3662  kern(ix1,iy2) += 1.0;
3663  sum += 1.0;
3664  }
3665  if ((ix2 < m_num_x) && (iy2 < m_num_y)) {
3666  kern(ix2,iy2) += 1.0;
3667  sum += 1.0;
3668  }
3669  }
3670  }
3671  }
3672 
3673  } // endif: handled disk kernel
3674 
3675  // Handle Gaussian kernel
3676  else if (gammalib::toupper(kernel) == "GAUSSIAN") {
3677 
3678  // Compute the exponential normalisation
3679  double norm = -0.5 / (par * par);
3680 
3681  // Fill kernel
3682  for (int ix1 = 0, ix2 = m_num_x; ix1 < m_num_x; ++ix1, --ix2) {
3683  double x = ix1 * dx;
3684  double xqs = x * x;
3685  for (int iy1 = 0, iy2 = m_num_y; iy1 < m_num_y; ++iy1, --iy2) {
3686  double y = iy1 * dy;
3687  double value = std::exp(norm * (xqs + y*y));
3688  kern(ix1,iy1) += value;
3689  sum += value;
3690  if (ix2 < m_num_x) {
3691  kern(ix2,iy1) += value;
3692  sum += value;
3693  }
3694  if (iy2 < m_num_y) {
3695  kern(ix1,iy2) += value;
3696  sum += value;
3697  }
3698  if ((ix2 < m_num_x) && (iy2 < m_num_y)) {
3699  kern(ix2,iy2) += value;
3700  sum += value;
3701  }
3702  }
3703  }
3704 
3705  } // endif: handled Gaussian kernel
3706 
3707  // ... otherwise throw an exception
3708  else {
3709  std::string msg = "Invalid kernel type \""+kernel+"\" specified. "
3710  "Please specify one of \"DISK\", \"GAUSSIAN\".";
3712  }
3713 
3714  // Normalise kernel
3715  if (sum > 0.0) {
3716  for (int ix = 0; ix < m_num_x; ++ix) {
3717  for (int iy = 0; iy < m_num_y; ++iy) {
3718  kern(ix,iy) /= sum;
3719  }
3720  }
3721  }
3722 
3723  // Return kernel
3724  return kern;
3725 }
3726 
3727 
3728 /*==========================================================================
3729  = =
3730  = Friends =
3731  = =
3732  ==========================================================================*/
3733 
3734 /***********************************************************************//**
3735  * @brief Computes square root of sky map elements
3736  *
3737  * @param[in] map Sky map.
3738  * @return Sky map containing the square root of every element.
3739  ***************************************************************************/
3740 GSkyMap sqrt(const GSkyMap& map)
3741 {
3742  // Initialise result vector
3743  GSkyMap result(map);
3744 
3745  // Loop over all maps
3746  for (int i = 0; i < map.nmaps(); ++i) {
3747 
3748  // Loop over all bins
3749  for (int j = 0; j < map.npix(); ++j) {
3750 
3751  // Get the content from the bin
3752  double content = map(j,i);
3753 
3754  // Throw an exception if content is negative
3755  if (content < 0.0) {
3756  std::string msg = "Negative value encountered. "
3757  "Cannot take the sqrt from a negative "
3758  "value";
3759  throw GException::invalid_value(G_SQRT, msg);
3760  }
3761 
3762  // Set content of the result map
3763  result(j,i) = std::sqrt(content);
3764 
3765  } // endfor: looped over all bins
3766 
3767  } // endfor: looped over all maps
3768 
3769  // Return sky map
3770  return result;
3771 }
3772 
3773 
3774 /***********************************************************************//**
3775  * @brief Computes the natural logarithm of sky map elements
3776  *
3777  * @param[in] map Sky map.
3778  * @return Sky map containing the natural logarithm of every element.
3779  ***************************************************************************/
3780 GSkyMap log(const GSkyMap& map)
3781 {
3782  // Initialise result vector
3783  GSkyMap result(map);
3784 
3785  // Loop over all maps
3786  for (int i = 0; i < map.nmaps(); ++i) {
3787 
3788  // Loop over all bins
3789  for (int j = 0; j < map.npix(); ++j) {
3790 
3791  // Get the content from the bin
3792  double content = map(j,i);
3793 
3794  // Throw an exception if content is negative
3795  if (content <= 0.0) {
3796  std::string msg = "Negative or zero value encountered. "
3797  "Cannot calculate the logarithm of a "
3798  "negative or zero value";
3799  throw GException::invalid_value(G_LOG, msg);
3800  }
3801 
3802  // Set content of the result map
3803  result(j,i) = std::log(content);
3804 
3805  } // endfor: looped over all bins
3806 
3807  } // endfor: looped over all maps
3808 
3809  // Return sky map
3810  return result;
3811 }
3812 
3813 
3814 /***********************************************************************//**
3815  * @brief Computes the base 10 logarithm of sky map elements
3816  *
3817  * @param[in] map Sky map.
3818  * @return Sky map containing the base 10 logarithm of every element.
3819  ***************************************************************************/
3820 GSkyMap log10(const GSkyMap& map)
3821 {
3822  // Initialise result vector
3823  GSkyMap result(map);
3824 
3825  // Loop over all maps
3826  for (int i = 0; i < map.nmaps(); ++i) {
3827 
3828  // Loop over all bins
3829  for (int j = 0; j < map.npix(); ++j) {
3830 
3831  // Get the content from the bin
3832  double content = map(j,i);
3833 
3834  // Throw an exception if content is negative
3835  if (content <= 0.0) {
3836  std::string msg = "Negative or zero value encountered. "
3837  "Cannot calculate the logarithm of a "
3838  "negative or zero value";
3839  throw GException::invalid_value(G_LOG10, msg);
3840  }
3841 
3842  // Set content of the result map
3843  result(j,i) = std::log10(content);
3844 
3845  } // endfor: looped over all bins
3846 
3847  } // endfor: looped over all maps
3848 
3849  // Return sky map
3850  return result;
3851 }
3852 
3853 
3854 /***********************************************************************//**
3855  * @brief Computes the absolute value of sky map elements
3856  *
3857  * @param[in] map Sky map.
3858  * @return Sky map containing the absolute value of every element.
3859  ***************************************************************************/
3860 GSkyMap abs(const GSkyMap& map)
3861 {
3862  // Initialise result vector
3863  GSkyMap result(map);
3864 
3865  // Loop over all maps
3866  for (int i = 0; i < map.nmaps(); ++i) {
3867 
3868  // Loop over all bins
3869  for (int j = 0; j < map.npix(); ++j) {
3870 
3871  // Get the content from the bin
3872  double content = map(j,i);
3873 
3874  // Set content of the result map
3875  result(j,i) = std::abs(content);
3876 
3877  } // endfor: looped over all bins
3878 
3879  } // endfor: looped over all maps
3880 
3881  // Return sky map
3882  return result;
3883 }
3884 
3885 
3886 /***********************************************************************//**
3887  * @brief Computes the sign value of sky map elements
3888  *
3889  * @param[in] map Sky map.
3890  * @return Sky map containing the sign value of every pixel.
3891  *
3892  * This method returns a sky map filled with a value of 1 if the pixel
3893  * is positive, a value of -1 if the pixel is negative or a value of 0
3894  * if the pixel is 0.
3895  ***************************************************************************/
3896 GSkyMap sign(const GSkyMap& map)
3897 {
3898  // Initialise result vector
3899  GSkyMap result(map);
3900 
3901  // Loop over all maps
3902  for (int i = 0; i < map.nmaps(); ++i) {
3903 
3904  // Loop over all bins
3905  for (int j = 0; j < map.npix(); ++j) {
3906 
3907  // Get the content from the bin
3908  double content = map(j,i);
3909 
3910  // Handle the 3 cases
3911  if (content < 0.0) {
3912  result(j,i) = -1.0;
3913  }
3914  else if (content > 0.0) {
3915  result(j,i) = +1.0;
3916  }
3917  else {
3918  result(j,i) = 0.0;
3919  }
3920 
3921  } // endfor: looped over all bins
3922 
3923  } // endfor: looed over all maps
3924 
3925  // Return sky map
3926  return result;
3927 }
3928 
3929 
3930 /***********************************************************************//**
3931  * @brief Clips map at given value
3932  *
3933  * @param[in] map Sky map.
3934  * @param[in] thresh Threshold value.
3935  * @return Sky map containing the value in the original map if >= thresh,
3936  * otherwise thresh
3937  ***************************************************************************/
3938 GSkyMap clip(const GSkyMap& map, const double& thresh)
3939 {
3940  // Initialise result vector
3941  GSkyMap result(map);
3942 
3943  // Loop over all maps
3944  for (int i = 0; i < map.nmaps(); ++i) {
3945 
3946  // Loop over all bins
3947  for (int j = 0; j < map.npix(); ++j) {
3948 
3949  // Get the content from the bin
3950  double content = map(j,i);
3951 
3952  // If content is below threshold, set to thresh
3953  if (content < thresh) {
3954  result(j,i) = thresh;
3955  }
3956 
3957  // ... otherwise keep content
3958  else {
3959  result(j,i) = content;
3960  }
3961 
3962  } // endfor: looped over all bins
3963 
3964  } // endfor: looped over all maps
3965 
3966  // Return sky map
3967  return result;
3968 }
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:3554
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
double dist_deg(const GSkyDir &dir) const
Compute angular distance between sky directions in degrees.
Definition: GSkyDir.hpp:284
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:1057
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:367
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:828
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:415
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:218
int m_num_pixels
Number of pixels (used for pixel allocation)
Definition: GSkyMap.hpp:217
void read_healpix(const GFitsTable &table)
Read Healpix data from FITS table.
Definition: GSkyMap.cpp:2890
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:854
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1170
std::vector< int > m_shape
Shape of the maps.
Definition: GSkyMap.hpp:221
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:219
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:2836
const double onethird
Definition: GMath.hpp:50
#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:222
double sum(const GVector &vector)
Computes vector sum.
Definition: GVector.cpp:908
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:882
int m_num_y
Number of pixels in y direction (only 2D)
Definition: GSkyMap.hpp:220
#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:2550
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:391
#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:1008
Abstract FITS image base class definition.
void publish(const std::string &name="") const
Publish sky map.
Definition: GSkyMap.cpp:2636
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:3503
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:2673
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:2770
#define G_OP_UNARY_ADD
Definition: GSkyMap.cpp:59
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1275
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:905
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:379
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:1296
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:3211
GFitsHDU * write(GFits &file, const std::string &extname="") const
Write sky map into FITS file.
Definition: GSkyMap.cpp:2583
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:223
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:928
int & index1(void)
Access index 1.
Definition: GBilinear.hpp:102
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1191
bool m_contained
Sky direction is contained in map.
Definition: GSkyMap.hpp:227
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:3110
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
GSkyRegionCircle region_circle(void) const
Return sky region circle that encloses the sky map.
Definition: GSkyMap.cpp:2287
void free_members(void)
Delete class members.
Definition: GSkyMap.cpp:2797
#define G_DIR2INX
Definition: GSkyMap.cpp:69
bool m_hascache
Cache is valid.
Definition: GSkyMap.hpp:226
Abstract world coordinate system base class.
Definition: GWcs.hpp:51
double & weight4(void)
Access weight 4.
Definition: GBilinear.hpp:186
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition: GFilename.hpp:189
const GSkyDir & centre(void) const
Return circular region centre.
void clear(void)
Clear sky direction.
Definition: GSkyDir.cpp:164
GSkyDir m_last_dir
Last sky direction.
Definition: GSkyMap.hpp:228
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:229
GFitsBinTable * create_healpix_hdu(void) const
Create FITS HDU containing Healpix data.
Definition: GSkyMap.cpp:3150
virtual double real(const int &row, const int &inx=0) const =0
int ndim(void) const
Returns dimension of maps.
Definition: GSkyMap.hpp:403
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:3529
void read_wcs(const GFitsImage &image)
Read WCS image from FITS HDU.
Definition: GSkyMap.cpp:3019
FITS binary table class.
bool overlaps_circle(const GSkyRegionCircle &region) const
Checks whether a circular region overlaps with this map.
Definition: GSkyMap.cpp:3440
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:3938
#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:269
std::string toupper(const std::string &s)
Convert string to upper case.
Definition: GTools.cpp:831
#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:351
#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:1149
void save(const GFilename &filename, const bool &clobber=false) const
Save sky map into FITS file.
Definition: GSkyMap.cpp:2497
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:1033
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:2397
const int & npix(void) const
Returns number of pixels.
Definition: GSkyMap.hpp:335
#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:2743
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:1065
#define G_SHAPE
Definition: GSkyMap.cpp:58
GNdarray smooth_kernel(const std::string &kernel, const double &par) const
Return smoothing kernel.
Definition: GSkyMap.cpp:3624
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:1212
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:415
FITS table abstract base class interface definition.