From 84cf03b48002cc4d2c43f753456cd63bcec11c88 Mon Sep 17 00:00:00 2001 From: Robert Butora <robert.butora@inaf.it> Date: Sat, 10 May 2025 14:59:14 +0200 Subject: [PATCH] cutout: makes POS optional (volib updated 0.9.8: NONE enum removed instead parsers return null if not found in input) --- .../engine/src/common/src/ast_frameset.cpp | 965 +++++++++--------- data-access/servlet/pom.xml | 2 +- .../servlet/src/main/java/common/Coord.java | 2 +- .../main/java/common/json/JsonEncoder.java | 1 - 4 files changed, 486 insertions(+), 484 deletions(-) diff --git a/data-access/engine/src/common/src/ast_frameset.cpp b/data-access/engine/src/common/src/ast_frameset.cpp index 2996470..1ea2271 100644 --- a/data-access/engine/src/common/src/ast_frameset.cpp +++ b/data-access/engine/src/common/src/ast_frameset.cpp @@ -964,132 +964,135 @@ overlap_ranges ast::frameset::overlap(coordinates coord) log_bounds("XX BOUNDS header: ", naxes, low, up); - /* overwrite bounds for sky-axis and spec axis if given by coord input */ - - int ix_lon = m_sky_lon_axis - 1; //astGetI(m_hdr_fs,"LonAxis") - 1; - int ix_lat = m_sky_lat_axis - 1; //astGetI(m_hdr_fs,"LatAxis") - 1; - //if ( !astOK ) - // throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetI( Lon/LatAxis )")); - - LOG_STREAM << "lon/lat axes indeces (zero-based): " - << to_string(ix_lon) << " " << to_string(ix_lat) << endl; - - set_bounds_to_rect(coord); - - low[ix_lon] = D2R * (coord.lon_deg - coord.dlon_deg/2.0); - low[ix_lat] = D2R * (coord.lat_deg - coord.dlat_deg/2.0); - up[ix_lon] = D2R * (coord.lon_deg + coord.dlon_deg/2.0); - 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)) - { - /* FIXME assumes m_frm_fs is VELO and unit is km/s How to align - if set(specsystem not called) ? */ - - string unit_key{"Unit("+to_string(m_spec_axis)+")"}; - const char * cunit = astGetC(m_hdr_fs, unit_key.c_str()); - string unit(cunit); - - string std_of_rest_key{"StdOfRest("+to_string(m_spec_axis)+")"}; - const char * cstd_of_rest = astGetC(m_hdr_fs, std_of_rest_key.c_str()); - string std_of_rest(cstd_of_rest); - - LOG_STREAM << "SpectralAxisUnit: " << unit << endl; - LOG_STREAM << "SpectralAxis StdOfRest: " << std_of_rest << endl; - - /* FIXME ast-overlap computation breaks if bands do not overlap - suspect AST bug ? */ - if((( up[m_spec_axis-1] < coord.vl_kmps) && ( up[m_spec_axis-1] < coord.vu_kmps)) || - ((low[m_spec_axis-1] > coord.vl_kmps) && (low[m_spec_axis-1] > coord.vu_kmps))) - { - /* set values only to get correct debug print for client */ - low[m_spec_axis-1] = coord.vl_kmps; - up [m_spec_axis-1] = coord.vu_kmps; - log_bounds("XX BOUNDS client: ", naxes, low, up); - LOG_STREAM << "XX BOUNDS no overlap in spectrum axis, returning ov_code=1" << endl; - return overlap_ranges{1, pix_ranges }; - } - else // at least partial overlap -> cut coord to min/max of header values - { - /* FIXME if() is needed becuase if coord bounds bigger then header's -> overlap yields 1 - suspect bug in AST? */ - if(low[m_spec_axis-1] < coord.vl_kmps) low[m_spec_axis-1] = coord.vl_kmps; - if(up [m_spec_axis-1] > coord.vu_kmps) up [m_spec_axis-1] = coord.vu_kmps; - } - } - - if(m_has_time_axis && (coord.timesys != timesystem::NONE)) - { - // FIXME see if() in spec bounds above: test verify and eventually removea comments - /*if(low[m_time_axis-1] < coord.time_value[0])*/ low[m_time_axis-1] = coord.time_value[0]; - /*if(up [m_time_axis-1] > coord.time_value[1])*/ up [m_time_axis-1] = coord.time_value[1]; - } - - if(m_has_stokes_axis && (!coord.pol.empty())) - { - // FIXME implement properly ranges 1,...4 and -1,...-8 (correspond to ranges I...V and RR...YY respectively) - low[m_stokes_axis-1] = min_pol_state(coord.pol); - up [m_stokes_axis-1] = max_pol_state(coord.pol); - } - - log_bounds("XX BOUNDS client: ", naxes, low, up); - - /* FIXME ignored coord.shape; add circle later : metters for ovelap-code (but not for cut, which is always rect) */ - AstRegion * crgn = (AstRegion *)astBox( wcsbox, 1, low, up , (AstRegion*)AST__NULL," "); - if ( !astOK ) - throw runtime_error(failed_with_status(__FILE__,__LINE__,"astBox( coord )")); - - astInvert(pix2wcs); - AstCmpRegion * crgn_pix = (AstCmpRegion*)astMapRegion( crgn, pix2wcs, pixfrm ); - if ( !astOK ) - throw runtime_error(failed_with_status(__FILE__,__LINE__,"astMapRegion( client-region WCS -> PIX by header-mapping )")); - - /* DBG only */ - double cllow[ AXES_CNT ]; - double clup[ AXES_CNT ]; - astGetRegionBounds(crgn_pix, cllow, clup); - if ( !astOK ) - throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetRegionBounds")); - - log_bounds("XX BOUNDS client_pix: ", naxes, cllow, clup); - /* DBG only */ - - - /* overlap */ - - /* FIXME returns "5"= exact macth however swapping coords returns seemingly correct overlap code ; bug ? - * FIXME TODO workround: if 5 returned -> swap regions and use that value - * accept 5 only if returned twice: original oreder and swapped both return 5 */ - //int ov_code = astOverlap( crgn_pix, pixbox ); - int ov_code = astOverlap( pixbox, crgn_pix ); - if ( !astOK ) - throw runtime_error(failed_with_status(__FILE__,__LINE__,"astOverlap( header, caller-coord )")); - - /* if overlap, calc pixel ranges */ - - - bool no_overlap = ((ov_code == 1) || (ov_code == 6)); - if(!no_overlap) - { - AstCmpRegion * pixOverlap = astCmpRegion( pixbox, crgn_pix, AST__AND, " "); - if ( !astOK || (pixOverlap == AST__NULL) ) - throw runtime_error(failed_with_status(__FILE__,__LINE__,"astCmpRegion( header AND coord )")); - - double lbnd[AXES_CNT]; - double ubnd[AXES_CNT]; - astGetRegionBounds(pixOverlap, lbnd, ubnd); - if ( !astOK ) - throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetRegionBounds( pixOverlap )")); - - log_bounds("XX BOUNDS overlap_pix: ", m_NAXIS, lbnd, ubnd); - - int ix; - for(ix=0; ix<m_NAXIS; ix++) - { - pix_ranges.push_back({lbnd[ix],ubnd[ix],m_axis_type[ix]}); - } - } - - return overlap_ranges{ov_code, pix_ranges }; + /* overwrite bounds for sky- spec- time- pol-axis if given by coord */ + + if(m_has_skyframe && (coord.skysys != skysystem::NONE)) + { + int ix_lon = m_sky_lon_axis - 1; //astGetI(m_hdr_fs,"LonAxis") - 1; + int ix_lat = m_sky_lat_axis - 1; //astGetI(m_hdr_fs,"LatAxis") - 1; + //if ( !astOK ) + // throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetI( Lon/LatAxis )")); + + LOG_STREAM << "lon/lat axes indeces (zero-based): " + << to_string(ix_lon) << " " << to_string(ix_lat) << endl; + + set_bounds_to_rect(coord); + + low[ix_lon] = D2R * (coord.lon_deg - coord.dlon_deg/2.0); + low[ix_lat] = D2R * (coord.lat_deg - coord.dlat_deg/2.0); + up[ix_lon] = D2R * (coord.lon_deg + coord.dlon_deg/2.0); + 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)) + { + /* FIXME assumes m_frm_fs is VELO and unit is km/s How to align - if set(specsystem not called) ? */ + + string unit_key{"Unit("+to_string(m_spec_axis)+")"}; + const char * cunit = astGetC(m_hdr_fs, unit_key.c_str()); + string unit(cunit); + + string std_of_rest_key{"StdOfRest("+to_string(m_spec_axis)+")"}; + const char * cstd_of_rest = astGetC(m_hdr_fs, std_of_rest_key.c_str()); + string std_of_rest(cstd_of_rest); + + LOG_STREAM << "SpectralAxisUnit: " << unit << endl; + LOG_STREAM << "SpectralAxis StdOfRest: " << std_of_rest << endl; + + /* FIXME ast-overlap computation breaks if bands do not overlap - suspect AST bug ? */ + if((( up[m_spec_axis-1] < coord.vl_kmps) && ( up[m_spec_axis-1] < coord.vu_kmps)) || + ((low[m_spec_axis-1] > coord.vl_kmps) && (low[m_spec_axis-1] > coord.vu_kmps))) + { + /* set values only to get correct debug print for client */ + low[m_spec_axis-1] = coord.vl_kmps; + up [m_spec_axis-1] = coord.vu_kmps; + log_bounds("XX BOUNDS client: ", naxes, low, up); + LOG_STREAM << "XX BOUNDS no overlap in spectrum axis, returning ov_code=1" << endl; + return overlap_ranges{1, pix_ranges }; + } + else // at least partial overlap -> cut coord to min/max of header values + { + /* FIXME if() is needed becuase if coord bounds bigger then header's -> overlap yields 1 - suspect bug in AST? */ + if(low[m_spec_axis-1] < coord.vl_kmps) low[m_spec_axis-1] = coord.vl_kmps; + if(up [m_spec_axis-1] > coord.vu_kmps) up [m_spec_axis-1] = coord.vu_kmps; + } + } + + if(m_has_time_axis && (coord.timesys != timesystem::NONE)) + { + // FIXME see if() in spec bounds above: test verify and eventually removea comments + /*if(low[m_time_axis-1] < coord.time_value[0])*/ low[m_time_axis-1] = coord.time_value[0]; + /*if(up [m_time_axis-1] > coord.time_value[1])*/ up [m_time_axis-1] = coord.time_value[1]; + } + + if(m_has_stokes_axis && (!coord.pol.empty())) + { + // FIXME implement properly ranges 1,...4 and -1,...-8 (correspond to ranges I...V and RR...YY respectively) + low[m_stokes_axis-1] = min_pol_state(coord.pol); + up [m_stokes_axis-1] = max_pol_state(coord.pol); + } + + log_bounds("XX BOUNDS client: ", naxes, low, up); + + /* FIXME ignored coord.shape; add circle later : metters for ovelap-code (but not for cut, which is always rect) */ + AstRegion * crgn = (AstRegion *)astBox( wcsbox, 1, low, up , (AstRegion*)AST__NULL," "); + if ( !astOK ) + throw runtime_error(failed_with_status(__FILE__,__LINE__,"astBox( coord )")); + + astInvert(pix2wcs); + AstCmpRegion * crgn_pix = (AstCmpRegion*)astMapRegion( crgn, pix2wcs, pixfrm ); + if ( !astOK ) + throw runtime_error(failed_with_status(__FILE__,__LINE__,"astMapRegion( client-region WCS -> PIX by header-mapping )")); + + /* DBG only */ + double cllow[ AXES_CNT ]; + double clup[ AXES_CNT ]; + astGetRegionBounds(crgn_pix, cllow, clup); + if ( !astOK ) + throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetRegionBounds")); + + log_bounds("XX BOUNDS client_pix: ", naxes, cllow, clup); + /* DBG only */ + + + /* overlap */ + + /* FIXME returns "5"= exact macth however swapping coords returns seemingly correct overlap code ; bug ? + * FIXME TODO workround: if 5 returned -> swap regions and use that value + * accept 5 only if returned twice: original oreder and swapped both return 5 */ + //int ov_code = astOverlap( crgn_pix, pixbox ); + int ov_code = astOverlap( pixbox, crgn_pix ); + if ( !astOK ) + throw runtime_error(failed_with_status(__FILE__,__LINE__,"astOverlap( header, caller-coord )")); + + /* if overlap, calc pixel ranges */ + + + bool no_overlap = ((ov_code == 1) || (ov_code == 6)); + if(!no_overlap) + { + AstCmpRegion * pixOverlap = astCmpRegion( pixbox, crgn_pix, AST__AND, " "); + if ( !astOK || (pixOverlap == AST__NULL) ) + throw runtime_error(failed_with_status(__FILE__,__LINE__,"astCmpRegion( header AND coord )")); + + double lbnd[AXES_CNT]; + double ubnd[AXES_CNT]; + astGetRegionBounds(pixOverlap, lbnd, ubnd); + if ( !astOK ) + throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetRegionBounds( pixOverlap )")); + + log_bounds("XX BOUNDS overlap_pix: ", m_NAXIS, lbnd, ubnd); + + int ix; + for(ix=0; ix<m_NAXIS; ix++) + { + pix_ranges.push_back({lbnd[ix],ubnd[ix],m_axis_type[ix]}); + } + } + + return overlap_ranges{ov_code, pix_ranges }; } #endif @@ -1098,193 +1101,193 @@ 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 - // and multiplied by the desired precision in ULPs (units in the last place) - return std::fabs(x-y) <= std::numeric_limits<double>::epsilon() * std::fabs(x+y) * ulp - // unless the result is subnormal - - || std::fabs(x-y) < std::numeric_limits<double>::min(); - - /* also note std::numeric_limits<double>::epsilon() cannot be used with a value having a magnitude different of 1.0. - * Instead, std::nextafter can be used to get the epsilon of a value with any magnitude. - * double nextafter (double x , double y ); - * returns the next representable value after x in the direction of y. - */ + // the machine epsilon has to be scaled to the magnitude of the values used + // and multiplied by the desired precision in ULPs (units in the last place) + return std::fabs(x-y) <= std::numeric_limits<double>::epsilon() * std::fabs(x+y) * ulp + // unless the result is subnormal + + || std::fabs(x-y) < std::numeric_limits<double>::min(); + + /* also note std::numeric_limits<double>::epsilon() cannot be used with a value having a magnitude different of 1.0. + * Instead, std::nextafter can be used to get the epsilon of a value with any magnitude. + * double nextafter (double x , double y ); + * returns the next representable value after x in the direction of y. + */ } /* 2**N */ int my_pow_int (int x, int p) { - int i = 1; - for (int j = 1; j <= p; j++) i *= x; - return i; + int i = 1; + for (int j = 1; j <= p; j++) i *= x; + return i; } std::vector<point2d> ast::frameset::sky_vertices(void) { - LOG_trace(__func__); + LOG_trace(__func__); - const int naxes = astGetI(m_hdr_fs, "Naxes"); - if ( !astOK ) - throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetI( Naxes )")); + const int naxes = astGetI(m_hdr_fs, "Naxes"); + if ( !astOK ) + throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetI( Naxes )")); - my_assert(naxes >= m_NAXIS, __FILE__,__LINE__,"AST-naxes is less then FITS-NAXIS"); + my_assert(naxes >= m_NAXIS, __FILE__,__LINE__,"AST-naxes is less then FITS-NAXIS"); - const int ncoord_in = naxes; - const int ncoord_out = naxes; // FIXME for now we assume Nin = Nout + const int ncoord_in = naxes; + const int ncoord_out = naxes; // FIXME for now we assume Nin = Nout - /* generate vertices in form (0.5, NAXISi+0.5) in 'naxes' dimensions */ + /* generate vertices in form (0.5, NAXISi+0.5) in 'naxes' dimensions */ - /* NOTE 1: lon lat coordinates can be identified only in WCS-coordinate due to - * possible permutations and transformations (rotation), so in pixel-grid - * we need to create complete vertex-array for each dimension */ + /* NOTE 1: lon lat coordinates can be identified only in WCS-coordinate due to + * possible permutations and transformations (rotation), so in pixel-grid + * we need to create complete vertex-array for each dimension */ - /* NOTE 2: this generates unordered set of vertices. - * If we could generate set of vertexes ordered by polar angle (from center of the 4-vertices) - * then polygons could be generated without later re-ordering. - * Is it possible at least in specific cases ? */ + /* NOTE 2: this generates unordered set of vertices. + * If we could generate set of vertexes ordered by polar angle (from center of the 4-vertices) + * then polygons could be generated without later re-ordering. + * Is it possible at least in specific cases ? */ - const int npoint = my_pow_int(2,m_NAXIS); + const int npoint = my_pow_int(2,m_NAXIS); - LOG_STREAM << "VERT Nin Nout: " << ncoord_in << " " << ncoord_out << " npoint: " << npoint << endl; + LOG_STREAM << "VERT Nin Nout: " << ncoord_in << " " << ncoord_out << " npoint: " << npoint << endl; - double vecs_in [ncoord_in ][npoint]; - double vecs_out[ncoord_out][npoint]; + double vecs_in [ncoord_in ][npoint]; + double vecs_out[ncoord_out][npoint]; - const double MIN_VAL = 0.5; + const double MIN_VAL = 0.5; - int ic,ip; - for(ip=0; ip < npoint; ip++) - { - for(ic=0; ic < ncoord_in; ic++) - { - const int bit = 1 << ic; + int ic,ip; + for(ip=0; ip < npoint; ip++) + { + for(ic=0; ic < ncoord_in; ic++) + { + const int bit = 1 << ic; - if(ic < m_NAXIS) - vecs_in[ic][ip] = (bit & ip) ? MIN_VAL : (MIN_VAL + m_NAXISn[ic]); - else - vecs_in[ic][ip] = (bit & ip) ? MIN_VAL : (MIN_VAL + 1); - } - } - - for(ip=0; ip<npoint; ip++) - { - LOG_STREAM << "Pi[" << ip << "]: "; - for(ic=0; ic<ncoord_in; ic++) - LOG_STREAM << " " << vecs_in[ic][ip]; - LOG_STREAM << endl; - } + if(ic < m_NAXIS) + vecs_in[ic][ip] = (bit & ip) ? MIN_VAL : (MIN_VAL + m_NAXISn[ic]); + else + vecs_in[ic][ip] = (bit & ip) ? MIN_VAL : (MIN_VAL + 1); + } + } + for(ip=0; ip<npoint; ip++) + { + LOG_STREAM << "Pi[" << ip << "]: "; + for(ic=0; ic<ncoord_in; ic++) + LOG_STREAM << " " << vecs_in[ic][ip]; + LOG_STREAM << endl; + } - /* transform the generated vertices in pixels to WCS and identify lon, lat axis */ - double * ptr_in [ncoord_in ]; - double * ptr_out[ncoord_out]; + /* transform the generated vertices in pixels to WCS and identify lon, lat axis */ - int ix; - for(ix=0; ix<ncoord_in ; ix++) ptr_in [ix] = &(vecs_in [ix][0]); - for(ix=0; ix<ncoord_out; ix++) ptr_out[ix] = &(vecs_out[ix][0]); + double * ptr_in [ncoord_in ]; + double * ptr_out[ncoord_out]; + int ix; + for(ix=0; ix<ncoord_in ; ix++) ptr_in [ix] = &(vecs_in [ix][0]); + for(ix=0; ix<ncoord_out; ix++) ptr_out[ix] = &(vecs_out[ix][0]); - astTranP(m_hdr_fs, npoint, ncoord_in, (const double**)ptr_in, 1, ncoord_out, ptr_out); - if ( !astOK ) - throw runtime_error(failed_with_status(__FILE__,__LINE__,"astTranP(pix -> wcs vertices)")); + astTranP(m_hdr_fs, npoint, ncoord_in, (const double**)ptr_in, 1, ncoord_out, ptr_out); + if ( !astOK ) + throw runtime_error(failed_with_status(__FILE__,__LINE__,"astTranP(pix -> wcs vertices)")); - const int ixlon = m_sky_lon_axis - 1; //astGetI(m_hdr_fs,"LonAxis") - 1; - const int ixlat = m_sky_lat_axis - 1; //astGetI(m_hdr_fs,"LatAxis") - 1; - //if ( !astOK ) - // throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetI( Lon/LatAxis )")); - LOG_STREAM << "VERT ix_lon ix_lat: " << ixlon << " " << ixlat << endl; + const int ixlon = m_sky_lon_axis - 1; //astGetI(m_hdr_fs,"LonAxis") - 1; + const int ixlat = m_sky_lat_axis - 1; //astGetI(m_hdr_fs,"LatAxis") - 1; + //if ( !astOK ) + // throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetI( Lon/LatAxis )")); + LOG_STREAM << "VERT ix_lon ix_lat: " << ixlon << " " << ixlat << endl; - /* select the 4 sky-vertices in case Naxes > 2 */ - double vertlon[npoint]; - double vertlat[npoint]; + /* select the 4 sky-vertices in case Naxes > 2 */ - int in = 0; - if(npoint <= 4) - { - for(ip=0; ip<npoint; ip++) - { - vertlon[ip] = ptr_out[ixlon][ip]; - vertlat[ip] = ptr_out[ixlat][ip]; - in++; - } - } - else - { - vertlon[in] = ptr_out[ixlon][0]; - vertlat[in] = ptr_out[ixlat][0]; - in++; + double vertlon[npoint]; + double vertlat[npoint]; - LOG_STREAM << "Po[0] cnt: <ref>" << endl; + int in = 0; + if(npoint <= 4) + { + for(ip=0; ip<npoint; ip++) + { + vertlon[ip] = ptr_out[ixlon][ip]; + vertlat[ip] = ptr_out[ixlat][ip]; + in++; + } + } + else + { + vertlon[in] = ptr_out[ixlon][0]; + vertlat[in] = ptr_out[ixlat][0]; + in++; - /* ip=0 used as reference */ - for(ip=1; ip<npoint; ip++) - { - int ic; - int cnt = 0; - for(ic=0; ic<ncoord_out; ic++) - { - if( (ic != ixlon) && (ic != ixlat) && (ptr_out[ic][ip] == ptr_out[ic][0]) ) cnt++; - /* asking for exact bit-pattern match on floating point numbers ptr_out[][] possible - * because axes are independent and so the same machine-code computation - * is performed on them which must yield the same result down to last bit */ - } + LOG_STREAM << "Po[0] cnt: <ref>" << endl; - bool coord_match = (cnt == (ncoord_out - 2)); + /* ip=0 used as reference */ + for(ip=1; ip<npoint; ip++) + { + int ic; + int cnt = 0; + for(ic=0; ic<ncoord_out; ic++) + { + if( (ic != ixlon) && (ic != ixlat) && (ptr_out[ic][ip] == ptr_out[ic][0]) ) cnt++; + /* asking for exact bit-pattern match on floating point numbers ptr_out[][] possible + * because axes are independent and so the same machine-code computation + * is performed on them which must yield the same result down to last bit */ + } - LOG_STREAM << "Po[" << ip << "] cnt: " << cnt << endl; + bool coord_match = (cnt == (ncoord_out - 2)); + + LOG_STREAM << "Po[" << ip << "] cnt: " << cnt << endl; - if(coord_match) - { - vertlon[in] = ptr_out[ixlon][ip]; - vertlat[in] = ptr_out[ixlat][ip]; - in++; - } + if(coord_match) + { + vertlon[in] = ptr_out[ixlon][ip]; + vertlat[in] = ptr_out[ixlat][ip]; + in++; + } #if 0 - double lon = ptr_out[ixlon][ip]; - double lat = ptr_out[ixlat][ip]; - - int ii; - bool found = false; - for(ii=0; ii<in; ii++) - { - int ulp = 2; // precision in units in the last place - if( almost_equal(vertlon[ii], lon, ulp) && almost_equal(vertlat[ii], lat, ulp) ) - { - found = true; - break; - } - } - - if(!found) - { - vertlon[in] = ptr_out[ixlon][ip]; - vertlat[in] = ptr_out[ixlat][ip]; - in++; - } + double lon = ptr_out[ixlon][ip]; + double lat = ptr_out[ixlat][ip]; + + int ii; + bool found = false; + for(ii=0; ii<in; ii++) + { + int ulp = 2; // precision in units in the last place + if( almost_equal(vertlon[ii], lon, ulp) && almost_equal(vertlat[ii], lat, ulp) ) + { + found = true; + break; + } + } + + if(!found) + { + vertlon[in] = ptr_out[ixlon][ip]; + vertlat[in] = ptr_out[ixlat][ip]; + in++; + } #endif - } - } + } + } - my_assert((in==4), __FILE__,__LINE__,"expected 4 vertices, but found " + to_string(in) + " from npoint: " + to_string(npoint)); + my_assert((in==4), __FILE__,__LINE__,"expected 4 vertices, but found " + to_string(in) + " from npoint: " + to_string(npoint)); - vector<point2d> ret_vert; + vector<point2d> ret_vert; - for(ix=0; ix<in; ix++) - { - ret_vert.push_back(point2d{(R2D)*vertlon[ix], (R2D)*vertlat[ix]}); + for(ix=0; ix<in; ix++) + { + ret_vert.push_back(point2d{(R2D)*vertlon[ix], (R2D)*vertlat[ix]}); - LOG_STREAM << "VERTx[" << ix << "] (deg): " << (R2D)*vertlon[ix] << " " << (R2D)*vertlat[ix] << endl; - } + LOG_STREAM << "VERTx[" << ix << "] (deg): " << (R2D)*vertlon[ix] << " " << (R2D)*vertlat[ix] << endl; + } - return ret_vert; + return ret_vert; } @@ -1294,204 +1297,204 @@ std::vector<point2d> ast::frameset::sky_vertices(void) void ast::frameset::assert_valid_state(void) { - LOG_trace(__func__); - - /* check NAXIS */ - - my_assert((m_NAXIS >= 0)&&(m_NAXIS < AXES_CNT), __FILE__,__LINE__, ": NAXIS must be > 0 and <= " + to_string(AXES_CNT)); - int ix; - for(ix=0; ix<m_NAXIS;ix++) - my_assert((m_NAXISn[ix] >= 0), __FILE__,__LINE__, ": NAXIS["+to_string(ix)+"]: " + to_string(m_NAXISn[ix])); - - /* check spec-axis identification */ - - if(m_has_specframe) - my_assert( (m_spec_axis >= 1)/*&&(m_spec_axis <=m_NAXIS) FIXME should be Naxes of m_hdr_fs */, - __FILE__,__LINE__, ": m_spec_axis is " + to_string(m_spec_axis) - + " and must be 1 .. " + to_string(m_NAXIS) ); - - if(m_has_skyframe) - { - my_assert( (m_sky_lon_axis >= 1), - __FILE__,__LINE__, ": m_sky_lon_axis is " + to_string(m_sky_lon_axis) - + " and must be 1 .. " + to_string(m_NAXIS) ); - my_assert( (m_sky_lat_axis >= 1), - __FILE__,__LINE__, ": m_sky_lat_axis is " + to_string(m_sky_lat_axis) - + " and must be 1 .. " + to_string(m_NAXIS) ); - } - - /* check AstFrameSet */ - - my_assert((m_hdr_fs != AST__NULL), __FILE__,__LINE__, ": AstFrameSet is NULL"); - /* FIXME add check internal to AstFrameSet: - * - FrameSet must have at least 2 Frames - * - AST__BASE Frame must have only GRID domain <-- verify this - * - AST__CURRENT Frame must have SKY domain - * - if AST__CURRENT Frame has also SpecFrame -> its axis must match m_spec_axis - * - if spectral present and is VELO unit should be km/s FIXME see in overlap - */ + LOG_trace(__func__); + + /* check NAXIS */ + + my_assert((m_NAXIS >= 0)&&(m_NAXIS < AXES_CNT), __FILE__,__LINE__, ": NAXIS must be > 0 and <= " + to_string(AXES_CNT)); + int ix; + for(ix=0; ix<m_NAXIS;ix++) + my_assert((m_NAXISn[ix] >= 0), __FILE__,__LINE__, ": NAXIS["+to_string(ix)+"]: " + to_string(m_NAXISn[ix])); + + /* check spec-axis identification */ + + if(m_has_specframe) + my_assert( (m_spec_axis >= 1)/*&&(m_spec_axis <=m_NAXIS) FIXME should be Naxes of m_hdr_fs */, + __FILE__,__LINE__, ": m_spec_axis is " + to_string(m_spec_axis) + + " and must be 1 .. " + to_string(m_NAXIS) ); + + if(m_has_skyframe) + { + my_assert( (m_sky_lon_axis >= 1), + __FILE__,__LINE__, ": m_sky_lon_axis is " + to_string(m_sky_lon_axis) + + " and must be 1 .. " + to_string(m_NAXIS) ); + my_assert( (m_sky_lat_axis >= 1), + __FILE__,__LINE__, ": m_sky_lat_axis is " + to_string(m_sky_lat_axis) + + " and must be 1 .. " + to_string(m_NAXIS) ); + } + + /* check AstFrameSet */ + + my_assert((m_hdr_fs != AST__NULL), __FILE__,__LINE__, ": AstFrameSet is NULL"); + /* FIXME add check internal to AstFrameSet: + * - FrameSet must have at least 2 Frames + * - AST__BASE Frame must have only GRID domain <-- verify this + * - AST__CURRENT Frame must have SKY domain + * - if AST__CURRENT Frame has also SpecFrame -> its axis must match m_spec_axis + * - if spectral present and is VELO unit should be km/s FIXME see in overlap + */ } std::ostream& operator<<(std::ostream& out, const ast::frameset& p) { - return p.serialize(out); + return p.serialize(out); } std::ostream& ast::frameset::serialize(std::ostream& ostrm) const { - ostrm << "ast::frameset: "; - - ostrm << "NAXIS["; - int ix; - for(ix=0; ix<(m_NAXIS-1); ix++) ostrm << m_NAXISn[ix] << " "; - ostrm << m_NAXISn[m_NAXIS-1]; - ostrm << "]"; - - ostrm << " Ast:"; - - const int n_frame = astGetI(m_hdr_fs,"Nframe"); - if ( !astOK ) - throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetI(Nframe)")); - - const char * domain = astGetC(m_hdr_fs,"Domain"); - if ( !astOK ) - throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetC(Domain)")); - - ostrm << "Domain(" << domain << ") "; - - - for(ix=0; ix<n_frame; ix++) - { - AstFrame * fr = (AstFrame*)astGetFrame(m_hdr_fs, ix+1); - - const char *astclass = astGetC( fr, "Class" ); - //const char *title = astGetC( fr, "Title" ); - const char *domain = astGetC( fr, "Domain" ); - //const char *system = astGetC( fr, "System" ); - const int naxes = astGetI(fr,"Naxes"); - if ( !astOK ) - throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetX()")); - - - ostrm << " "; - ostrm << astclass<< "[" << ix+1 << "]{"; - ostrm << "Domain(" << domain << ")"; - - ostrm << " Axes["; - int axis; - for(axis=1; axis<=naxes; axis++) - { - const string domain_key{"Domain(" + to_string(axis) + ")"}; - const string symbol_key{"Symbol(" + to_string(axis) + ")"}; - const string system_key{"System(" + to_string(axis) + ")"}; - const string unit_key{"Unit(" + to_string(axis) + ")"}; - - const char * domain = astGetC(fr, domain_key.c_str()); - const char * symbol = astGetC(fr, symbol_key.c_str()); - const char * system = astGetC(fr, system_key.c_str()); - const char * unit = astGetC(fr, unit_key.c_str()); - if ( !astOK ) - throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetX()")); - - ostrm <<domain << " " << symbol << " " << system << " " << unit; - if(axis != naxes) ostrm << " | "; - } - ostrm << "]"; - ostrm << "}"; - } - - - ostrm << " has specframe: " << m_has_specframe; - if(m_has_specframe) - { - ostrm << " at axis (one-based) " << to_string(m_spec_axis) << " "; - - string unit_key{"Unit("+to_string(m_spec_axis)+")"}; - const char * cunit = astGetC(m_hdr_fs, unit_key.c_str()); - string unit(cunit); - LOG_STREAM << "SpecUnit(" << unit << ") "; - - string std_of_rest_key{"StdOfRest("+to_string(m_spec_axis)+")"}; - const char * cstd_of_rest = astGetC(m_hdr_fs, std_of_rest_key.c_str()); - string std_of_rest(cstd_of_rest); - LOG_STREAM << "StdOfRest(" << std_of_rest << ")"; - } - - ostrm << " | "; - - ostrm << "has stokes axis: " << m_has_stokes_axis; - if(m_has_stokes_axis) - { - ostrm << " at axis (one-based) " << to_string(m_stokes_axis) << " "; - } - - ostrm << " | "; - - ostrm << "has time axis: " << m_has_time_axis; - if(m_has_time_axis) - { - ostrm << " at axis (one-based) " << to_string(m_time_axis) << " "; - - string system_key{"System("+to_string(m_time_axis)+")"}; - const char * csystem = astGetC(m_hdr_fs, system_key.c_str()); - string system(csystem); - LOG_STREAM << "System(" << system << ")"; - - string unit_key{"Unit("+to_string(m_time_axis)+")"}; - const char * cunit = astGetC(m_hdr_fs, unit_key.c_str()); - string unit(cunit); - LOG_STREAM << "TimeUnit(" << unit << ") "; - } - - ostrm << " | "; - - ostrm << "has skyframe: " << m_has_skyframe; - if(m_has_skyframe) - { - ostrm << " [lon,lat] at axes (one-based) [" - << to_string(m_sky_lon_axis) << "," - << to_string(m_sky_lat_axis) << "] "; - - string system_key_lon{"System("+to_string(m_sky_lon_axis)+")"}; - const char * csystem_lon = astGetC(m_hdr_fs, system_key_lon.c_str()); - string system_lon(csystem_lon); - LOG_STREAM << "System[Lon](" << system_lon << ") "; - - string system_key_lat{"System("+to_string(m_sky_lat_axis)+")"}; - const char * csystem_lat = astGetC(m_hdr_fs, system_key_lat.c_str()); - string system_lat(csystem_lat); - LOG_STREAM << "System[Lat](" << system_lat << ") "; - - string unit_key_lon{"Unit("+to_string(m_sky_lon_axis)+")"}; - const char * cunit_lon = astGetC(m_hdr_fs, unit_key_lon.c_str()); - string unit_lon(cunit_lon); - LOG_STREAM << "Unit[Lon](" << unit_lon << ") "; - - string unit_key_lat{"Unit("+to_string(m_sky_lat_axis)+")"}; - const char * cunit_lat = astGetC(m_hdr_fs, unit_key_lat.c_str()); - string unit_lat(cunit_lat); - LOG_STREAM << "Unit[Lat](" << unit_lat << ") "; - } - - ostrm << " axes types: >"; - //for(char cc : m_axis_type) ostrm << cc; - for(int ix=0; ix<AXES_CNT;ix++) ostrm << m_axis_type[ix]; - ostrm << "<"; - - ostrm << endl; - - return ostrm; + ostrm << "ast::frameset: "; + + ostrm << "NAXIS["; + int ix; + for(ix=0; ix<(m_NAXIS-1); ix++) ostrm << m_NAXISn[ix] << " "; + ostrm << m_NAXISn[m_NAXIS-1]; + ostrm << "]"; + + ostrm << " Ast:"; + + const int n_frame = astGetI(m_hdr_fs,"Nframe"); + if ( !astOK ) + throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetI(Nframe)")); + + const char * domain = astGetC(m_hdr_fs,"Domain"); + if ( !astOK ) + throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetC(Domain)")); + + ostrm << "Domain(" << domain << ") "; + + + for(ix=0; ix<n_frame; ix++) + { + AstFrame * fr = (AstFrame*)astGetFrame(m_hdr_fs, ix+1); + + const char *astclass = astGetC( fr, "Class" ); + //const char *title = astGetC( fr, "Title" ); + const char *domain = astGetC( fr, "Domain" ); + //const char *system = astGetC( fr, "System" ); + const int naxes = astGetI(fr,"Naxes"); + if ( !astOK ) + throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetX()")); + + + ostrm << " "; + ostrm << astclass<< "[" << ix+1 << "]{"; + ostrm << "Domain(" << domain << ")"; + + ostrm << " Axes["; + int axis; + for(axis=1; axis<=naxes; axis++) + { + const string domain_key{"Domain(" + to_string(axis) + ")"}; + const string symbol_key{"Symbol(" + to_string(axis) + ")"}; + const string system_key{"System(" + to_string(axis) + ")"}; + const string unit_key{"Unit(" + to_string(axis) + ")"}; + + const char * domain = astGetC(fr, domain_key.c_str()); + const char * symbol = astGetC(fr, symbol_key.c_str()); + const char * system = astGetC(fr, system_key.c_str()); + const char * unit = astGetC(fr, unit_key.c_str()); + if ( !astOK ) + throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetX()")); + + ostrm <<domain << " " << symbol << " " << system << " " << unit; + if(axis != naxes) ostrm << " | "; + } + ostrm << "]"; + ostrm << "}"; + } + + + ostrm << " has specframe: " << m_has_specframe; + if(m_has_specframe) + { + ostrm << " at axis (one-based) " << to_string(m_spec_axis) << " "; + + string unit_key{"Unit("+to_string(m_spec_axis)+")"}; + const char * cunit = astGetC(m_hdr_fs, unit_key.c_str()); + string unit(cunit); + LOG_STREAM << "SpecUnit(" << unit << ") "; + + string std_of_rest_key{"StdOfRest("+to_string(m_spec_axis)+")"}; + const char * cstd_of_rest = astGetC(m_hdr_fs, std_of_rest_key.c_str()); + string std_of_rest(cstd_of_rest); + LOG_STREAM << "StdOfRest(" << std_of_rest << ")"; + } + + ostrm << " | "; + + ostrm << "has stokes axis: " << m_has_stokes_axis; + if(m_has_stokes_axis) + { + ostrm << " at axis (one-based) " << to_string(m_stokes_axis) << " "; + } + + ostrm << " | "; + + ostrm << "has time axis: " << m_has_time_axis; + if(m_has_time_axis) + { + ostrm << " at axis (one-based) " << to_string(m_time_axis) << " "; + + string system_key{"System("+to_string(m_time_axis)+")"}; + const char * csystem = astGetC(m_hdr_fs, system_key.c_str()); + string system(csystem); + LOG_STREAM << "System(" << system << ")"; + + string unit_key{"Unit("+to_string(m_time_axis)+")"}; + const char * cunit = astGetC(m_hdr_fs, unit_key.c_str()); + string unit(cunit); + LOG_STREAM << "TimeUnit(" << unit << ") "; + } + + ostrm << " | "; + + ostrm << "has skyframe: " << m_has_skyframe; + if(m_has_skyframe) + { + ostrm << " [lon,lat] at axes (one-based) [" + << to_string(m_sky_lon_axis) << "," + << to_string(m_sky_lat_axis) << "] "; + + string system_key_lon{"System("+to_string(m_sky_lon_axis)+")"}; + const char * csystem_lon = astGetC(m_hdr_fs, system_key_lon.c_str()); + string system_lon(csystem_lon); + LOG_STREAM << "System[Lon](" << system_lon << ") "; + + string system_key_lat{"System("+to_string(m_sky_lat_axis)+")"}; + const char * csystem_lat = astGetC(m_hdr_fs, system_key_lat.c_str()); + string system_lat(csystem_lat); + LOG_STREAM << "System[Lat](" << system_lat << ") "; + + string unit_key_lon{"Unit("+to_string(m_sky_lon_axis)+")"}; + const char * cunit_lon = astGetC(m_hdr_fs, unit_key_lon.c_str()); + string unit_lon(cunit_lon); + LOG_STREAM << "Unit[Lon](" << unit_lon << ") "; + + string unit_key_lat{"Unit("+to_string(m_sky_lat_axis)+")"}; + const char * cunit_lat = astGetC(m_hdr_fs, unit_key_lat.c_str()); + string unit_lat(cunit_lat); + LOG_STREAM << "Unit[Lat](" << unit_lat << ") "; + } + + ostrm << " axes types: >"; + //for(char cc : m_axis_type) ostrm << cc; + for(int ix=0; ix<AXES_CNT;ix++) ostrm << m_axis_type[ix]; + ostrm << "<"; + + ostrm << endl; + + return ostrm; } void ast::frameset::log_NAXISn(void) { - LOG_STREAM << "NAXISn[AXES_CNT]: "; - int ix; - for(ix=0; ix<AXES_CNT;ix++) LOG_STREAM << m_NAXISn[ix] << " "; - LOG_STREAM << endl; + LOG_STREAM << "NAXISn[AXES_CNT]: "; + int ix; + for(ix=0; ix<AXES_CNT;ix++) LOG_STREAM << m_NAXISn[ix] << " "; + LOG_STREAM << endl; } @@ -1499,42 +1502,42 @@ void ast::frameset::log_NAXISn(void) void ast::frameset::log_warnings(const AstFitsChan * fchan) { - LOG_trace(__func__); + LOG_trace(__func__); - // reports warnings of the last astRead or astWrite invocation + // reports warnings of the last astRead or astWrite invocation - AstKeyMap *warnings; + AstKeyMap *warnings; #define KEYLEN (15) - char key[ KEYLEN ]; - short int iwarn; - const char *message; - - if( astGetI( fchan, "ReportLevel" ) > 0 ) - { - warnings = (AstKeyMap*)astWarnings( fchan ); - - if( warnings && astOK ) - { - LOG_STREAM << "The following warnings were issued" << endl; - - iwarn = 1; - while( astOK ) - { - snprintf( key,KEYLEN, "Warning_%d", iwarn++ ); - key[KEYLEN-1]='\0'; - if( astMapGet0C( warnings, key, &message ) ) - { - LOG_STREAM << string{key} << " - " << string{message} << endl; - } - else - { - break; - } - } - } - /* FIXME why clear ? */ - astClearStatus; - } + char key[ KEYLEN ]; + short int iwarn; + const char *message; + + if( astGetI( fchan, "ReportLevel" ) > 0 ) + { + warnings = (AstKeyMap*)astWarnings( fchan ); + + if( warnings && astOK ) + { + LOG_STREAM << "The following warnings were issued" << endl; + + iwarn = 1; + while( astOK ) + { + snprintf( key,KEYLEN, "Warning_%d", iwarn++ ); + key[KEYLEN-1]='\0'; + if( astMapGet0C( warnings, key, &message ) ) + { + LOG_STREAM << string{key} << " - " << string{message} << endl; + } + else + { + break; + } + } + } + /* FIXME why clear ? */ + astClearStatus; + } } diff --git a/data-access/servlet/pom.xml b/data-access/servlet/pom.xml index 94280f4..12c56cd 100644 --- a/data-access/servlet/pom.xml +++ b/data-access/servlet/pom.xml @@ -65,7 +65,7 @@ <dependency> <groupId>vo</groupId> <artifactId>vlkb-volib</artifactId> - <version>0.9.7</version> + <version>0.9.8</version> </dependency> <!-- dependency> diff --git a/data-access/servlet/src/main/java/common/Coord.java b/data-access/servlet/src/main/java/common/Coord.java index 7f364e2..6f74dac 100644 --- a/data-access/servlet/src/main/java/common/Coord.java +++ b/data-access/servlet/src/main/java/common/Coord.java @@ -7,7 +7,7 @@ class Coord String skySystem; // FIXME enum ICRS | GALACTIC String shape = "CIRCLE"; // FIXME enum CIRCLE | RECT | POLYGON FIXME replace RECT -> RANGE - String specSystem; // FIXME enum VELO_LSRK | WAVE_Barycentric | NONE + String specSystem; // FIXME enum VELO_LSRK | WAVE_Barycentric Pos pos; Band band; diff --git a/data-access/servlet/src/main/java/common/json/JsonEncoder.java b/data-access/servlet/src/main/java/common/json/JsonEncoder.java index d0cffff..99d6ed8 100644 --- a/data-access/servlet/src/main/java/common/json/JsonEncoder.java +++ b/data-access/servlet/src/main/java/common/json/JsonEncoder.java @@ -27,7 +27,6 @@ public class JsonEncoder // NOTE: implementation of NULL: do not put into json - // (alternatively could put JSON.NULL as "pos" value or use pos.system.NONE) public void add(Pos pos) { -- GitLab