Justine Devin / gammalib

Commit 07fd4b9edbc4a756dc719caa4938abb0b89b5e97

Authored by Jürgen Knödlseder
1 parent 6e9d8a89

Handle wrap around in X axis and add check of skymap size after reading. 

The wrap around in the X axis is now correctly dealt with in the interpolation operator.

Upon reading of the WCS map the size of the map is checked. Exceptions are thrown in case that the map is larger then 360 in X and larger than 180 in Y.
Inline Side-by-side
Showing 1 changed file with 66 additions and 27 deletions
src/sky/GSkymap.cpp
  View file @07fd4b9
... ... @@ -344,7 +344,7 @@ GSkymap& GSkymap::operator+=(const GSkymap& map)
344 344 std::string msg = "Mismatch of number of maps in skymap object"
345 345 " ("+gammalib::str(nmaps())+" maps in destination"
346 346 " map, "+gammalib::str(map.nmaps())+" in source"
347   - " map.";
  347 + " map).";
348 348 throw GException::invalid_value(G_OP_UNARY_ADD, msg);
349 349 }
350 350  
... ... @@ -418,7 +418,7 @@ GSkymap& GSkymap::operator-=(const GSkymap& map)
418 418 std::string msg = "Mismatch of number of maps in skymap object"
419 419 " ("+gammalib::str(nmaps())+" maps in destination"
420 420 " map, "+gammalib::str(map.nmaps())+" in source"
421   - " map.";
  421 + " map).";
422 422 throw GException::invalid_value(G_OP_UNARY_SUB, msg);
423 423 }
424 424  
... ... @@ -492,7 +492,7 @@ GSkymap& GSkymap::operator*=(const GSkymap& map)
492 492 std::string msg = "Mismatch of number of maps in skymap object"
493 493 " ("+gammalib::str(nmaps())+" maps in destination"
494 494 " map, "+gammalib::str(map.nmaps())+" in source"
495   - " map.";
  495 + " map).";
496 496 throw GException::invalid_value(G_OP_UNARY_MUL, msg);
497 497 }
498 498  
... ... @@ -567,7 +567,7 @@ GSkymap& GSkymap::operator/=(const GSkymap& map)
567 567 std::string msg = "Mismatch of number of maps in skymap object"
568 568 " ("+gammalib::str(nmaps())+" maps in destination"
569 569 " map, "+gammalib::str(map.nmaps())+" in source"
570   - " map.";
  570 + " map).";
571 571 throw GException::invalid_value(G_OP_UNARY_DIV, msg);
572 572 }
573 573  
... ... @@ -847,35 +847,53 @@ double GSkymap::operator()(const GSkyDir& dir, const int& map) const
847 847 // Continue only if pixel is within the map
848 848 if (m_contained) {
849 849  
850   - // Set left indices for interpolation. The left index is
851   - // comprised between 0 and npixels-2. By definition, the
852   - // right index is then the left index + 1
853   - int inx_x = int(pixel.x());
854   - int inx_y = int(pixel.y());
855   - if (inx_x < 0) {
856   - inx_x = 0;
  850 + // Signal that we have a wrap around in the x axis
  851 + double x_size = std::abs(m_num_x * static_cast<GWcs*>(m_proj)->cdelt(0));
  852 + bool x_wrap = (x_size > 359.99);
  853 +
  854 + // Set left and right indices for interpolation. The left
  855 + // index needs to non-negative and the right index needs
  856 + // to be not larger that the number of pixels. We treat
  857 + // here also the special case of wrap around in the x
  858 + // axis that may occur if we have an allsky map.
  859 + int inx_x_left = int(pixel.x());
  860 + int inx_y_left = int(pixel.y());
  861 + int inx_x_right = inx_x_left + 1;
  862 + int inx_y_right = inx_y_left + 1;
  863 + if (inx_x_left < 0) {
  864 + if (x_wrap) {
  865 + inx_x_left += m_num_x;
  866 + }
  867 + else {
  868 + inx_x_left = 0;
  869 + }
857 870 }
858   - else if (inx_x > m_num_x-2) {
859   - inx_x = m_num_x - 2;
  871 + if (inx_y_right >= m_num_x) {
  872 + if (x_wrap) {
  873 + inx_y_right -= m_num_x;
  874 + }
  875 + else {
  876 + inx_y_right = m_num_x - 1;
  877 + }
860 878 }
861   - if (inx_y < 0) {
862   - inx_y = 0;
  879 + if (inx_y_left < 0) {
  880 + inx_y_left = 0;
863 881 }
864   - else if (inx_y > m_num_y-2) {
865   - inx_y = m_num_y - 2;
  882 + if (inx_y_right >= m_num_y) {
  883 + inx_y_right = m_num_y - 1;
866 884 }
867 885  
868 886 // Set weighting factors for interpolation
869   - double wgt_x_right = (pixel.x() - inx_x);
  887 + double wgt_x_right = (pixel.x() - int(pixel.x()));
870 888 double wgt_x_left = 1.0 - wgt_x_right;
871   - double wgt_y_right = (pixel.y() - inx_y);
  889 + double wgt_y_right = (pixel.y() - int(pixel.y()));
872 890 double wgt_y_left = 1.0 - wgt_y_right;
873 891  
874 892 // Compute skymap pixel indices for bi-linear interpolation
875   - m_interpol.index1() = inx_x + inx_y * m_num_x;
876   - m_interpol.index2() = m_interpol.index1() + m_num_x;
877   - m_interpol.index3() = m_interpol.index1() + 1;
878   - m_interpol.index4() = m_interpol.index2() + 1;
  893 + m_interpol.index1() = inx_x_left + inx_y_left * m_num_x;
  894 + m_interpol.index2() = inx_x_left + inx_y_right * m_num_x;
  895 + m_interpol.index3() = inx_x_right + inx_y_left * m_num_x;
  896 + m_interpol.index4() = inx_x_right + inx_y_right * m_num_x;
879 897  
880 898 // Compute weighting factors for bi-linear interpolation
881 899 m_interpol.weight1() = wgt_x_left * wgt_y_left;
... ... @@ -1644,7 +1662,6 @@ void GSkymap::load(const std::string&amp; filename)
1644 1662  
1645 1663 // Load WCS map
1646 1664 read_wcs(static_cast<const GFitsImage&>(hdu));
1647   - //loaded = true;
1648 1665 break;
1649 1666  
1650 1667 } // endfor: looped over HDUs
... ... @@ -2129,8 +2146,9 @@ void GSkymap::read_healpix(const GFitsTable&amp; table)
2129 2146 *
2130 2147 * @param[in] image FITS image.
2131 2148 *
2132   - * @exception GException::skymap_bad_image_dim
2133   - * WCS image has invalid dimension (naxis=2 or 3).
  2149 + * @exception GException::invalid_value
  2150 + * WCS image has an invalid dimension or covers a too large
  2151 + * range.
2134 2152 ***************************************************************************/
2135 2153 void GSkymap::read_wcs(const GFitsImage& image)
2136 2154 {
... ... @@ -2152,7 +2170,10 @@ void GSkymap::read_wcs(const GFitsImage&amp; image)
2152 2170 m_num_maps = image.naxes(2);
2153 2171 }
2154 2172 else {
2155   - throw GException::skymap_bad_image_dim(G_READ_WCS, image.naxis());
  2173 + std::string msg = "Skymap has a dimension of "+
  2174 + gammalib::str(image.naxis())+" but only dimensions "
  2175 + "of 2 or 3 are supported.";
  2176 + throw GException::invalid_value(G_READ_WCS, msg);
2156 2177 }
2157 2178 #if defined(G_READ_WCS_DEBUG)
2158 2179 std::cout << "m_num_x=" << m_num_x << std::endl;
... ... @@ -2189,6 +2210,24 @@ void GSkymap::read_wcs(const GFitsImage&amp; image)
2189 2210 }
2190 2211 }
2191 2212  
  2213 + // Check if the map is too large
  2214 + double x_size = std::abs(m_num_x * static_cast<GWcs*>(m_proj)->cdelt(0));
  2215 + double y_size = std::abs(m_num_y * static_cast<GWcs*>(m_proj)->cdelt(1));
  2216 + if (x_size > 360.001) {
  2217 + std::string msg = "Skymap covers "+gammalib::str(x_size)+" degrees "
  2218 + "in the X axis which results in ambiguous "
  2219 + "coordinate transformations. Please provide a "
  2220 + "skymap that does not cover more than 360 degrees.";
  2221 + throw GException::invalid_value(G_READ_WCS, msg);
  2222 + }
  2223 + if (y_size > 180.001) {
  2224 + std::string msg = "Skymap covers "+gammalib::str(y_size)+" degrees "
  2225 + "in the Y axis which results in ambiguous "
  2226 + "coordinate transformations. Please provide a "
  2227 + "skymap that does not cover more than 180 degrees.";
  2228 + throw GException::invalid_value(G_READ_WCS, msg);
  2229 + }
  2230 +
2192 2231 // Return
2193 2232 return;
2194 2233 }
... ...