Commit 125ff464 authored by Robert Butora's avatar Robert Butora
Browse files

cutout: fixes double->long conversion for pixel-grid in overlap calculation on domain=GRID

parent 2f6018fd

data-access/engine/src/common/include/ast4vl.hpp

+3 −10
Original line number Diff line number Diff line
@@ -5,7 +5,7 @@ 
#include <string>

#include <vector>



#include "cutout.hpp" // coordinate needed

#include "cutout.hpp" // uint_bounds, coordinate needed





struct point2d {double lon; double lat;};
@@ -21,17 +21,10 @@ struct Bounds 
   int naxis;

};



struct double_xy

{

   double x;

   double y;

   /*unsigned*/ char type;

};



struct overlap_ranges

{

   int ov_code;

   std::vector<double_xy> pixel_ranges;

   std::vector<uint_bounds> pixel_ranges;

};
@@ -47,7 +40,7 @@ std::vector<point2d> calc_skyvertices(std::string header, std::string skysys); 


std::vector<Bounds> calc_bounds(std::string header, std::string skysys_str, std::string specsys_str);



std::vector<uint_bounds> calc_overlap(const std::string header, const coordinates coord, int& ov_code);

overlap_ranges calc_overlap(const std::string header, const coordinates coord);



int header_modif_coordsys(const std::string& skysys_str, const std::string& specsys_str,

      const std::string& pathname);

data-access/engine/src/common/src/ast4vl.cpp

+3 −61
Original line number Diff line number Diff line
@@ -38,7 +38,7 @@ std::ostream& operator<<( std::ostream &out, struct uint_bounds const& p) 
std::ostream& operator<<( std::ostream &out, overlap_ranges const& p)

{

   out << p.ov_code;

   for(double_xy r : p.pixel_ranges) out  <<" ("<< r.x << ", " << r.y  << ")";

   for(uint_bounds r : p.pixel_ranges) out  <<" ("<< r.pix1 << ", " << r.pix2  << ")";



   return out;

}
@@ -106,7 +106,7 @@ std::vector<Bounds> calc_bounds(std::string header, std::string skysys, std::str 
 */





std::vector<uint_bounds> calc_overlap(const std::string header, const coordinates coord, int& ov_code)

overlap_ranges calc_overlap(const std::string header, const coordinates coord)

{

   LOG_trace(__func__);
@@ -124,67 +124,9 @@ std::vector<uint_bounds> calc_overlap(const std::string header, const coordinate 


   overlap_ranges pix_ranges = frm_set.overlap(coord);



   ov_code = pix_ranges.ov_code;



   LOG_STREAM << "ov-code & pix ranges[double]: " << pix_ranges << endl;



   /* convert to uint */



   vector<uint_bounds> uint_bounds_vec;



   LOG_STREAM << "uint_bounds: " << uint_bounds_vec.size()

              << "  pix_ranges:: " << pix_ranges.pixel_ranges.size()

              << endl;





   LOG_STREAM << "pix ranges[uint]:";

   int ix = 0;

   for(double_xy dbl_range : pix_ranges.pixel_ranges)

   {

      if(dbl_range.x < 0)

         throw out_of_range(string{__FILE__} + ":" + to_string(__LINE__)

               + " pixel axis from overlap x is negative " + to_string(dbl_range.x));



      if(dbl_range.y < 0)

         throw out_of_range(string{__FILE__} + ":" + to_string(__LINE__)

               + " pixel axis from overlap y is negative " + to_string(dbl_range.y));



      // conversion double -> uint: result must be: 1 <= result <= NAXISn

      // because NAXISn start with 1 (FITS standard) which corresponds to ASTlib's GRID-domain

      // FitsChan uses GRID Domain for FITS-pixel coords

      // in GRID first element represents a pixel-range: [0.5..1.5) (0.5 in pixel_1 ;  1.5 in pixel_2)

      //

      // dbl-range is checked to be within 0.5 .. NAXIS+0.5 in ast_frameset::overlap()

      if(dbl_range.x <= dbl_range.y)

      {

         long lowb = lround(dbl_range.x);

         long upb  = lround(dbl_range.y);

			long NAXISi = frm_set.get_naxis(ix++);



         if(lowb == (NAXISi + 1)) lowb--;

         if(upb  == (NAXISi + 1)) upb--;



         uint_bounds ui_range{lowb, upb, dbl_range.type};

         uint_bounds_vec.push_back(ui_range);

         LOG_STREAM << " " << ui_range;

      }

      else

      {

         long lowb = lround(dbl_range.y);

         long upb  = lround(dbl_range.x);

			long NAXISi = frm_set.get_naxis(ix++);



         if(lowb == (NAXISi + 1)) lowb--;

         if(upb  == (NAXISi + 1)) upb--;



         uint_bounds ui_range{lowb, upb, dbl_range.type};

         uint_bounds_vec.push_back(ui_range);

         LOG_STREAM << "s" << ui_range;

      }

   }

   LOG_STREAM << endl;



   return uint_bounds_vec;

   return pix_ranges;

}

data-access/engine/src/common/src/ast_frameset.cpp

+14 −13
Original line number Diff line number Diff line
@@ -147,11 +147,6 @@ ast::frameset::frameset(string header) 
}





long ast::frameset::get_naxis(int ix)

{

   return m_NAXISn[ix];  

}





bool is_specdomain(AstFrame * frm)

{
@@ -882,8 +877,12 @@ void ast::frameset::set_bounds_to_rect(coordinates& coord) 
	}

}







long lround_grid(double num, long NAXISi)

{

   long lnum = lround(num);

   if(lnum == (NAXISi+1)) lnum--;

   return lnum;

}



/* calcs overlap in pixel coords */

overlap_ranges ast::frameset::overlap(coordinates coord)
@@ -994,7 +993,6 @@ overlap_ranges ast::frameset::overlap(coordinates coord) 
		up[ix_lat]  = D2R * (coord.lat_deg + coord.dlat_deg/2.0);

	}



	vector<double_xy> pix_ranges;



	if(m_has_specframe && (coord.specsys != specsystem::NONE))

	{
@@ -1020,7 +1018,9 @@ overlap_ranges ast::frameset::overlap(coordinates coord) 
			up [m_spec_axis-1] = coord.vu_kmps;

			log_bounds("BOUNDS client WCS: ", naxes, low, up);

			LOG_STREAM << "BOUNDS no overlap in spectrum axis, returning ov_code=1" << endl;

			return overlap_ranges{1, pix_ranges };



	      vector<uint_bounds> empty_pix_ranges;

			return overlap_ranges{1, empty_pix_ranges };

		}

		else // at least partial overlap -> cut coord to min/max of header values

		{
@@ -1079,7 +1079,7 @@ overlap_ranges ast::frameset::overlap(coordinates coord) 


	/* if overlap, calc pixel ranges */





	vector<uint_bounds> pix_ranges;

	bool no_overlap = ((ov_code == 1) || (ov_code == 6));

	if(!no_overlap)

	{
@@ -1118,7 +1118,9 @@ overlap_ranges ast::frameset::overlap(coordinates coord) 
                       + to_string(ix+1) + "] + 0.5 = " + to_string(p1[ix])

                       + ": " + to_string(ubnd[ix]));



			pix_ranges.push_back({lbnd[ix],ubnd[ix],m_axis_type[ix]});

			pix_ranges.push_back( { lround_grid(lbnd[ix],m_NAXISn[ix]),

                                 lround_grid(ubnd[ix],m_NAXISn[ix]),

                                 m_axis_type[ix]} );

		}

	}
@@ -1128,7 +1130,6 @@ overlap_ranges ast::frameset::overlap(coordinates coord) 








bool almost_equal(double x, double y, int ulp)

{

	// the machine epsilon has to be scaled to the magnitude of the values used

data-access/engine/src/common/src/ast_frameset.hpp

+1 −3
Original line number Diff line number Diff line
@@ -7,7 +7,7 @@ extern "C" { 
}



#include "cutout.hpp" // coordinates needed

#include "ast4vl.hpp" // struct decls needed

#include "ast4vl.hpp" // overlap_ranges,  struct decls needed



#include <array>

#include <vector>
@@ -24,8 +24,6 @@ class frameset 
      frameset(std::string header);

      ~frameset();



      long get_naxis(int ix);



      bool has_specframe(void);

      bool has_timeaxis(void);

data-access/engine/src/common/src/cutout.cpp

+3 −2
Original line number Diff line number Diff line
@@ -304,8 +304,9 @@ std::uintmax_t cutout_file( 
   }



   int ov_code;

   vector<uint_bounds> bounds = calc_overlap(header, coord, ov_code);

   //vector<uint_bounds> bounds = legacy::call_AST4VL_overlap(coord, header, ov_code);

   overlap_ranges ov_ranges = calc_overlap(header, coord);

   ov_code = ov_ranges.ov_code;

   vector<uint_bounds> bounds = ov_ranges.pixel_ranges;



   if((ov_code==1)||(ov_code==6))

      throw invalid_argument("given coordinates do not overlap with given fits-hdu area (ov_code = " + to_string(ov_code) +")");