Samuel Flis / gammalib

Commit 089f4f795dce95a6dd51dc85a9ef943e8d15a746

Authored by Jürgen Knödlseder
1 parent 23329634

Correct sky map smoothing 

The smoothing kernel covered before only two sectors, now all four sectors are covered. The smoothing now also works correctly in case that more than a single sky map is comprised in a GSkyMap object.
Inline Side-by-side
Showing 5 changed files with 32 additions and 9 deletions
ChangeLog
  View file @089f4f7
1   -2018-02-08
  1 +2018-02-13
2 2  
3 3 * Bug fix version 1.5.1 released
4 4 ================================
5 5  
  6 + Correct sky map smoothing.
  7 +
6 8  
7 9 2018-01-23
8 10  
... ...
1 1 New Features and Important Changes in GammaLib 1.5.1
2 2  
3   -8 February 2018
  3 +13 February 2018
4 4  
5 5  
6 6 1. Introduction
... ... @@ -83,7 +83,9 @@ None
83 83  
84 84 12. Sky module
85 85 --------------
86   -None
  86 +The smoothing kernel covered before only two sectors, now all four sectors
  87 +are covered. The smoothing now also works correctly in case that more than
  88 +a single sky map is comprised in a GSkyMap object.
87 89  
88 90  
89 91 13. Support module
... ...
doc/source/admin/release_history/1.5.1.rst
  View file @089f4f7
... ... @@ -10,7 +10,8 @@ GammaLib 1.5.1 is a bug fix release for GammaLib release 1.5.
10 10  
11 11 The following bugs were fixed:
12 12  
13   -* ...
  13 +* Correct sky map smoothing
  14 +
14 15  
15 16 The following improvements were made:
16 17  
... ...
src/sky/GSkyMap.cpp
  View file @089f4f7
... ... @@ -1866,11 +1866,11 @@ void GSkyMap::smooth(const std::string& kernel, const double& par)
1866 1866 GFft fft_kernel(smooth_kernel(kernel, par));
1867 1867  
1868 1868 // Loop over all sky maps
1869   - for (int i = 0; i < m_num_maps; ++i) {
  1869 + for (int k = 0; k < m_num_maps; ++k) {
1870 1870  
1871 1871 // Extract sky map
1872 1872 GNdarray array(m_num_x, m_num_y);
1873   - const double *src = m_pixels.data() + i*m_num_pixels;
  1873 + const double *src = m_pixels.data() + k*m_num_pixels;
1874 1874 double *dst = array.data();
1875 1875 for (int i = 0; i < m_num_pixels; ++i) {
1876 1876 *dst++ = *src++;
... ... @@ -1882,12 +1882,12 @@ void GSkyMap::smooth(const std::string&amp; kernel, const double&amp; par)
1882 1882 // Multiply FFT of sky map with FFT of kernel
1883 1883 GFft fft_smooth = fft_array * fft_kernel;
1884 1884  
1885   - // Background transform sky map
  1885 + // Backward transform sky map
1886 1886 GNdarray smooth = fft_smooth.backward();
1887 1887  
1888 1888 // Insert sky map
1889 1889 src = smooth.data();
1890   - dst = m_pixels.data() + i*m_num_pixels;
  1890 + dst = m_pixels.data() + k*m_num_pixels;
1891 1891 for (int i = 0; i < m_num_pixels; ++i) {
1892 1892 *dst++ = *src++;
1893 1893 }
... ... @@ -3187,7 +3187,7 @@ GNdarray GSkyMap::smooth_kernel(const std::string&amp; kernel,
3187 3187 const GWcs* wcs = dynamic_cast<const GWcs*>(m_proj);
3188 3188 if (wcs == NULL) {
3189 3189 std::string msg = "Sky map is not a WCS projection. Method is only "
3190   - "valid for WCS projects.";
  3190 + "valid for WCS projections.";
3191 3191 throw GException::invalid_argument(G_SMOOTH_KERNEL, msg);
3192 3192 }
3193 3193  
... ... @@ -3213,6 +3213,14 @@ GNdarray GSkyMap::smooth_kernel(const std::string&amp; kernel,
3213 3213 if (std::sqrt(xqs + y*y) <= par) {
3214 3214 kern(ix1,iy1) += 1.0;
3215 3215 sum += 1.0;
  3216 + if (ix2 < m_num_x) {
  3217 + kern(ix2,iy1) += 1.0;
  3218 + sum += 1.0;
  3219 + }
  3220 + if (iy2 < m_num_y) {
  3221 + kern(ix1,iy2) += 1.0;
  3222 + sum += 1.0;
  3223 + }
3216 3224 if ((ix2 < m_num_x) && (iy2 < m_num_y)) {
3217 3225 kern(ix2,iy2) += 1.0;
3218 3226 sum += 1.0;
... ... @@ -3238,6 +3246,14 @@ GNdarray GSkyMap::smooth_kernel(const std::string&amp; kernel,
3238 3246 double value = std::exp(norm * (xqs + y*y));
3239 3247 kern(ix1,iy1) += value;
3240 3248 sum += value;
  3249 + if (ix2 < m_num_x) {
  3250 + kern(ix2,iy1) += value;
  3251 + sum += value;
  3252 + }
  3253 + if (iy2 < m_num_y) {
  3254 + kern(ix1,iy2) += value;
  3255 + sum += value;
  3256 + }
3241 3257 if ((ix2 < m_num_x) && (iy2 < m_num_y)) {
3242 3258 kern(ix2,iy2) += value;
3243 3259 sum += value;
... ...
test/test_GSky.cpp
  View file @089f4f7
... ... @@ -1212,6 +1212,7 @@ void TestGSky::test_GSkyMap(void)
1212 1212 // Smooth map using DISK kernel
1213 1213 GSkyMap map_smooth1 = map_smooth;
1214 1214 map_smooth1.smooth("DISK", 0.2);
  1215 + map_smooth1.save("test_map_smooth_disk.fits", true);
1215 1216  
1216 1217 // Compute sum of smoothed sky map
1217 1218 double sum_smooth1 = 0.0;
... ... @@ -1223,6 +1224,7 @@ void TestGSky::test_GSkyMap(void)
1223 1224 // Smooth map using GAUSSIAN kernel
1224 1225 GSkyMap map_smooth2 = map_smooth;
1225 1226 map_smooth2.smooth("GAUSSIAN", 0.2);
  1227 + map_smooth2.save("test_map_smooth_gaussian.fits", true);
1226 1228  
1227 1229 // Compute sum of smoothed sky map
1228 1230 double sum_smooth2 = 0.0;
... ...