From e988d0e6af7e4daae39811409cefad60d2a34398 Mon Sep 17 00:00:00 2001
From: Robert Butora <robert.butora@inaf.it>
Date: Sat, 10 May 2025 16:47:11 +0200
Subject: [PATCH] cutout: improves debug-print in overlap
---
.../engine/src/common/include/cutout.hpp | 16 +-
.../engine/src/common/src/ast_frameset.cpp | 275 +++++++++---------
.../engine/src/common/src/ast_frameset.hpp | 1 +
.../engine/src/common/src/cutout_ostream.cpp | 12 +-
4 files changed, 155 insertions(+), 149 deletions(-)
diff --git a/data-access/engine/src/common/include/cutout.hpp b/data-access/engine/src/common/include/cutout.hpp
index 72adb2f..16670a5 100644
--- a/data-access/engine/src/common/include/cutout.hpp
+++ b/data-access/engine/src/common/include/cutout.hpp
@@ -68,14 +68,14 @@ struct coordinates
{
coordinates();
- skysystem skysys; // mandatory: "GALACTIC"|"ICRS"
- double lon_deg, lat_deg; // mandatory
- area shape; // mandatory: CIRCLE dlon=dlat=2xRadius, RECT dlon may differ from dlat
- double dlon_deg, dlat_deg; // mandatory
- std::vector<double> p_lon_deg; // mandatory polygon arrry
- std::vector<double> p_lat_deg; // mandatory polygon array
-
- specsystem specsys; // mandatory NONE VELO-LSRK WAVE-Bary
+ skysystem skysys; // mandatory: NONE GALACTIC ICRS
+ double lon_deg, lat_deg; // optional, check skysys == NONE
+ area shape; // optional: CIRCLE dlon=dlat=2xRadius, RECT dlon may differ dlat
+ double dlon_deg, dlat_deg; // optional
+ std::vector<double> p_lon_deg; // optional polygon arrry
+ std::vector<double> p_lat_deg; // optional polygon array
+
+ specsystem specsys; // mandatory NONE VELO_LSRK WAVE_Barycentric
double vl_kmps, vu_kmps; // optional, check specsystem == NONE
timesystem timesys;
diff --git a/data-access/engine/src/common/src/ast_frameset.cpp b/data-access/engine/src/common/src/ast_frameset.cpp
index 1ea2271..7731661 100644
--- a/data-access/engine/src/common/src/ast_frameset.cpp
+++ b/data-access/engine/src/common/src/ast_frameset.cpp
@@ -791,12 +791,17 @@ overlap_ranges ast::frameset::overlap_in_wcs(coordinates coord)
return overlap_ranges{ov_code, pix_ranges };
}
#else
-void log_bounds(string prefix, int length, double low[], double up[])
+void ast::frameset::log_bounds(string prefix, int length, double low[], double up[])
{
- int ix;
- LOG_STREAM << prefix << endl;
- for(ix=0; ix<length; ix++) LOG_STREAM << prefix << low[ix] << " " << up[ix]
- << " deg: " << R2D*low[ix] << " " << R2D*up[ix] << endl;
+ int ix;
+ double factor;
+ bool is_wcs = (prefix.find("WCS") != std::string::npos);
+ for(ix=0; ix<length; ix++)
+ {
+ factor = (is_wcs && (ix == (m_sky_lon_axis-1) || ix == (m_sky_lat_axis-1))) ? R2D : 1.0;
+ LOG_STREAM << prefix << factor*low[ix] << " " << factor*up[ix] << endl;
+ }
+ // << " deg: " << R2D*low[ix] << " " << R2D*up[ix] << endl;
}
@@ -804,72 +809,72 @@ void log_bounds(string prefix, int length, double low[], double up[])
void ast::frameset::set_bounds_to_rect(coordinates& coord)
{
- LOG_trace(__func__);
+ LOG_trace(__func__);
- switch(coord.shape)
- {
- case area::CIRCLE:
- case area::RECT:
- /* both set already */
- break;
- case area::POLYGON:
- {
- my_assert(coord.p_lon_deg.size()==coord.p_lat_deg.size(),
- __FILE__,__LINE__,"coord::p_lon and p_lat sizes differ");
-
- constexpr int int_max_value {std::numeric_limits<int>::max()};
- my_assert(coord.p_lon_deg.size() > int_max_value,
- __FILE__,__LINE__,"coord-polygon too long for astPolygon");
-
- int npnt = (int)coord.p_lon_deg.size();
- int dim = npnt;
- double points[2][dim];
- const double * pts = &(points[0][0]);
-
- LOG_STREAM << "polygon ";
- int ii;
- for(ii = 0; ii<dim; ii++)
- {
- points[0][ii] = D2R * coord.p_lon_deg[ii];
- points[1][ii] = D2R * coord.p_lat_deg[ii];
-
- LOG_STREAM << "(" << R2D * points[0][ii] << ", " << R2D * points[1][ii] << ")";
- }
- LOG_STREAM << endl;
-
- AstSkyFrame * sky_frm = (AstSkyFrame*)find_skyframe();
-
- my_assert(sky_frm != nullptr, __FILE__,__LINE__,"sky frame not found in header frameset");
-
- AstPolygon * astPoly = astPolygon(sky_frm, npnt, dim, pts, NULL, " ");
- if ( !astOK )
- throw runtime_error(failed_with_status(__FILE__,__LINE__,"astPolygon( npoints=" + to_string(npnt) + "...)"));
-
- double lbnd[2];
- double ubnd[2];
- astGetRegionBounds( (AstRegion*)astPoly, lbnd, ubnd );
- if ( !astOK )
- throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetRegionBounds(astPolygon ...)"));
- /* NOTE AST-manual:
- * lbnd < ubnd -> ok, note: if axis has no lower/upper limit largest-negative/largest-positive value is returned
- * lbnd = ubnd -> axis has one value
- * lbnd > ubnd -> region has no extent on that axis
- * lbnd = ubnd = AST__NULL -> bounds on that axis cannot be determined */
-
- LOG_STREAM << "polygon bounds lon: " << R2D * lbnd[0] << " .. " << R2D * ubnd[0] << endl;
- LOG_STREAM << "polygon bounds lat: " << R2D * lbnd[1] << " .. " << R2D * ubnd[1] << endl;
-
- // (mis)use RECT in coord to store bounds:
- coord.lon_deg = R2D * (lbnd[0] + ubnd[0])/2.0;
- coord.lat_deg = R2D * (lbnd[1] + ubnd[1])/2.0;
- coord.dlon_deg = R2D * (ubnd[0] - lbnd[0]);
- coord.dlat_deg = R2D * (ubnd[1] - lbnd[1]);
- }
- break;
+ switch(coord.shape)
+ {
+ case area::CIRCLE:
+ case area::RECT:
+ /* both set already */
+ break;
+ case area::POLYGON:
+ {
+ my_assert(coord.p_lon_deg.size()==coord.p_lat_deg.size(),
+ __FILE__,__LINE__,"coord::p_lon and p_lat sizes differ");
- default:
- my_assert(false, __FILE__,__LINE__,"coord::shape invalid");
- }
+ constexpr int int_max_value {std::numeric_limits<int>::max()};
+ my_assert(coord.p_lon_deg.size() > int_max_value,
+ __FILE__,__LINE__,"coord-polygon too long for astPolygon");
+
+ int npnt = (int)coord.p_lon_deg.size();
+ int dim = npnt;
+ double points[2][dim];
+ const double * pts = &(points[0][0]);
+
+ LOG_STREAM << "polygon ";
+ int ii;
+ for(ii = 0; ii<dim; ii++)
+ {
+ points[0][ii] = D2R * coord.p_lon_deg[ii];
+ points[1][ii] = D2R * coord.p_lat_deg[ii];
+
+ LOG_STREAM << "(" << R2D * points[0][ii] << ", " << R2D * points[1][ii] << ")";
+ }
+ LOG_STREAM << endl;
+
+ AstSkyFrame * sky_frm = (AstSkyFrame*)find_skyframe();
+
+ my_assert(sky_frm != nullptr, __FILE__,__LINE__,"sky frame not found in header frameset");
+
+ AstPolygon * astPoly = astPolygon(sky_frm, npnt, dim, pts, NULL, " ");
+ if ( !astOK )
+ throw runtime_error(failed_with_status(__FILE__,__LINE__,"astPolygon( npoints=" + to_string(npnt) + "...)"));
+
+ double lbnd[2];
+ double ubnd[2];
+ astGetRegionBounds( (AstRegion*)astPoly, lbnd, ubnd );
+ if ( !astOK )
+ throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetRegionBounds(astPolygon ...)"));
+ /* NOTE AST-manual:
+ * lbnd < ubnd -> ok, note: if axis has no lower/upper limit largest-negative/largest-positive value is returned
+ * lbnd = ubnd -> axis has one value
+ * lbnd > ubnd -> region has no extent on that axis
+ * lbnd = ubnd = AST__NULL -> bounds on that axis cannot be determined */
+
+ LOG_STREAM << "polygon bounds lon: " << R2D * lbnd[0] << " .. " << R2D * ubnd[0] << endl;
+ LOG_STREAM << "polygon bounds lat: " << R2D * lbnd[1] << " .. " << R2D * ubnd[1] << endl;
+
+ // (mis)use RECT in coord to store bounds:
+ coord.lon_deg = R2D * (lbnd[0] + ubnd[0])/2.0;
+ coord.lat_deg = R2D * (lbnd[1] + ubnd[1])/2.0;
+ coord.dlon_deg = R2D * (ubnd[0] - lbnd[0]);
+ coord.dlat_deg = R2D * (ubnd[1] - lbnd[1]);
+ }
+ break;
+
+ default:
+ my_assert(false, __FILE__,__LINE__,"coord::shape invalid");
+ }
}
@@ -878,95 +883,95 @@ void ast::frameset::set_bounds_to_rect(coordinates& coord)
/* calcs overlap in pixel coords */
overlap_ranges ast::frameset::overlap(coordinates coord)
{
- LOG_trace(__func__);
+ LOG_trace(__func__);
- /* ---------- get all header bounds ---------------------- */
+ /* ---------- get all header bounds ---------------------- */
- /* AstRegion * wcsbox = create_header_region();
- * NOTE: not used because it creates header-region based on NAXIS only, not Naxes:
- * FITS header defines number of axes as max of ( NAXIS,WCSAXIS and highest wcs-index in a wcs-key) */
+ /* AstRegion * wcsbox = create_header_region();
+ * NOTE: not used because it creates header-region based on NAXIS only, not Naxes:
+ * FITS header defines number of axes as max of ( NAXIS,WCSAXIS and highest wcs-index in a wcs-key) */
- AstFrame * wcsfrm = (AstFrame*)astGetFrame( m_hdr_fs, AST__CURRENT );
- if ( !astOK )
- throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetFrame( CURRENT )"));
+ AstFrame * wcsfrm = (AstFrame*)astGetFrame( m_hdr_fs, AST__CURRENT );
+ if ( !astOK )
+ throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetFrame( CURRENT )"));
- AstFrame * pixfrm = (AstFrame*)astGetFrame( m_hdr_fs, AST__BASE );
- if ( !astOK )
- throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetFrame( BASE )"));
+ AstFrame * pixfrm = (AstFrame*)astGetFrame( m_hdr_fs, AST__BASE );
+ if ( !astOK )
+ throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetFrame( BASE )"));
- AstMapping * pix2wcs = (AstMapping*)astGetMapping( m_hdr_fs, AST__BASE, AST__CURRENT );
- if ( !astOK )
- throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetMapping( BASE->CURRENT )"));
+ AstMapping * pix2wcs = (AstMapping*)astGetMapping( m_hdr_fs, AST__BASE, AST__CURRENT );
+ if ( !astOK )
+ throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetMapping( BASE->CURRENT )"));
- /* Create AstBox from NAXISn */
+ /* Create AstBox from NAXISn */
- double p1[ AXES_CNT ];
- double p2[ AXES_CNT ];
+ double p1[ AXES_CNT ];
+ double p2[ AXES_CNT ];
- my_assert(m_NAXIS <= AXES_CNT, __FILE__,__LINE__,
- "This build supports " + to_string(AXES_CNT) + " axes, but NAXIS is bigger : " + to_string(m_NAXIS));
+ my_assert(m_NAXIS <= AXES_CNT, __FILE__,__LINE__,
+ "This build supports " + to_string(AXES_CNT) + " axes, but NAXIS is bigger : " + to_string(m_NAXIS));
- int ix;
- for(ix = 0; ix < m_NAXIS; ix++)
- {
- p1[ix] = 0.5;
- p2[ix] = 0.5 + m_NAXISn[ix];
- }
+ int ix;
+ for(ix = 0; ix < m_NAXIS; ix++)
+ {
+ p1[ix] = 0.5;
+ p2[ix] = 0.5 + m_NAXISn[ix];
+ }
- int naxes = astGetI( m_hdr_fs, "Naxes" );
- if ( !astOK )
- throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetI( Naxes )"));
+ int naxes = astGetI( m_hdr_fs, "Naxes" );
+ if ( !astOK )
+ throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetI( Naxes )"));
- LOG_STREAM << "Naxes / NAXIS : " << naxes << " / " << m_NAXIS << endl;
+ LOG_STREAM << "Naxes / NAXIS : " << naxes << " / " << m_NAXIS << endl;
- my_assert(naxes <= AXES_CNT, __FILE__,__LINE__,
- "This build supports " + to_string(AXES_CNT)
- + " axes, but Naxes returned from astGetI(header-frameset, Naxes) is bigger : " + to_string(m_NAXIS));
+ my_assert(naxes <= AXES_CNT, __FILE__,__LINE__,
+ "This build supports " + to_string(AXES_CNT)
+ + " axes, but Naxes returned from astGetI(header-frameset, Naxes) is bigger : " + to_string(m_NAXIS));
- for(ix = m_NAXIS; ix < naxes; ix++)
- {
- p1[ix] = 0.5;
- p2[ix] = 1.5;
- }
+ for(ix = m_NAXIS; ix < naxes; ix++)
+ {
+ p1[ix] = 0.5;
+ p2[ix] = 1.5;
+ }
- AstBox * pixbox = astBox( pixfrm, 1, p1, p2, NULL, " " );
- if ( !astOK )
- throw runtime_error(failed_with_status(__FILE__,__LINE__,"astBox( NAXISn )"));
+ AstBox * pixbox = astBox( pixfrm, 1, p1, p2, NULL, " " );
+ if ( !astOK )
+ throw runtime_error(failed_with_status(__FILE__,__LINE__,"astBox( NAXISn )"));
- /* DEBUG only */
- double xlbnd[AXES_CNT];
- double xubnd[AXES_CNT];
- astGetRegionBounds(pixbox, xlbnd, xubnd);
- if ( !astOK )
- throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetRegionBounds( pixbox )"));
- log_bounds("XX BOUNDS header_pix: ", naxes, xlbnd, xubnd);
- /* DEBUG only */
+ /* DEBUG only */
+ double xlbnd[AXES_CNT];
+ double xubnd[AXES_CNT];
+ astGetRegionBounds(pixbox, xlbnd, xubnd);
+ if ( !astOK )
+ throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetRegionBounds( pixbox )"));
+ log_bounds("BOUNDS header PIX: ", naxes, xlbnd, xubnd);
+ /* DEBUG only */
- /* map header's pixels to WCS-region */
+ /* map header's pixels to WCS-region */
- AstRegion * wcsbox = (AstRegion*)astMapRegion( pixbox, pix2wcs, wcsfrm );
- if ( !astOK )
- throw runtime_error(failed_with_status(__FILE__,__LINE__,"astMapRegion( pixbox -> wcsbox )"));
+ AstRegion * wcsbox = (AstRegion*)astMapRegion( pixbox, pix2wcs, wcsfrm );
+ if ( !astOK )
+ throw runtime_error(failed_with_status(__FILE__,__LINE__,"astMapRegion( pixbox -> wcsbox )"));
- /* ------------ get bounds of the header region -------------------------- */
+ /* ------------ get bounds of the header region -------------------------- */
- double low[ AXES_CNT ];
- double up[ AXES_CNT ];
+ double low[ AXES_CNT ];
+ double up[ AXES_CNT ];
- astGetRegionBounds(wcsbox, low, up);
- if ( !astOK )
- throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetRegionBounds"));
+ astGetRegionBounds(wcsbox, low, up);
+ if ( !astOK )
+ throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetRegionBounds"));
- log_bounds("XX BOUNDS header: ", naxes, low, up);
+ log_bounds("BOUNDS header WCS: ", naxes, low, up);
- /* overwrite bounds for sky- spec- time- pol-axis if given by coord */
+ /* overwrite bounds for sky- spec- time- pol-axis if given by coord */
- if(m_has_skyframe && (coord.skysys != skysystem::NONE))
+ 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;
@@ -1008,8 +1013,8 @@ overlap_ranges ast::frameset::overlap(coordinates coord)
/* 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;
+ 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 };
}
else // at least partial overlap -> cut coord to min/max of header values
@@ -1034,7 +1039,7 @@ overlap_ranges ast::frameset::overlap(coordinates coord)
up [m_stokes_axis-1] = max_pol_state(coord.pol);
}
- log_bounds("XX BOUNDS client: ", naxes, low, up);
+ log_bounds("BOUNDS client WCS: ", 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," ");
@@ -1053,7 +1058,7 @@ overlap_ranges ast::frameset::overlap(coordinates coord)
if ( !astOK )
throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetRegionBounds"));
- log_bounds("XX BOUNDS client_pix: ", naxes, cllow, clup);
+ log_bounds("BOUNDS client PIX: ", naxes, cllow, clup);
/* DBG only */
@@ -1083,7 +1088,7 @@ overlap_ranges ast::frameset::overlap(coordinates coord)
if ( !astOK )
throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetRegionBounds( pixOverlap )"));
- log_bounds("XX BOUNDS overlap_pix: ", m_NAXIS, lbnd, ubnd);
+ log_bounds("BOUNDS overlap PIX: ", m_NAXIS, lbnd, ubnd);
int ix;
for(ix=0; ix<m_NAXIS; ix++)
diff --git a/data-access/engine/src/common/src/ast_frameset.hpp b/data-access/engine/src/common/src/ast_frameset.hpp
index b1c570a..73ec69b 100644
--- a/data-access/engine/src/common/src/ast_frameset.hpp
+++ b/data-access/engine/src/common/src/ast_frameset.hpp
@@ -45,6 +45,7 @@ class frameset
std::ostream& serialize(std::ostream& strm) const;
private:
+ void log_bounds(std::string prefix, int length, double low[], double up[]);
void log_warnings(const AstFitsChan * fchan);
AstRegion * create_header_region(void);
diff --git a/data-access/engine/src/common/src/cutout_ostream.cpp b/data-access/engine/src/common/src/cutout_ostream.cpp
index 3499df4..eb4a3d0 100644
--- a/data-access/engine/src/common/src/cutout_ostream.cpp
+++ b/data-access/engine/src/common/src/cutout_ostream.cpp
@@ -110,15 +110,15 @@ std::ostream& operator<<( std::ostream &out, struct coordinates const& p)
bool is_poly = (p.shape == area::POLYGON);
- out << to_string(p.skysys) << " : (" << p.lon_deg << ", " << p.lat_deg << ") "
- << to_string(p.shape) << " : ("
+ out << "POS["<< to_string(p.skysys) << "; (" << p.lon_deg << ", " << p.lat_deg << ") "
+ << to_string(p.shape) << ": ("
<< (is_poly ? poly_lon_str : to_string(p.dlon_deg)) << ", "
- << (is_poly ? poly_lat_str : to_string(p.dlat_deg)) << ") "
- << to_string(p.specsys) << " [km/s] : (" << p.vl_kmps << ", " << p.vu_kmps << ") "
+ << (is_poly ? poly_lat_str : to_string(p.dlat_deg)) << ")] "
+ << "BAND[" <<to_string(p.specsys) << "; (" << p.vl_kmps << ", " << p.vu_kmps << ")] "
// << "POS(" << p.pos << ") "
// << "BAND(" << to_string(p.bandsys) << "; " << p.band_value[0] << ", " << p.band_value[1] << ") "
- << "TIME(" << to_string(p.timesys) << "; " << p.time_value[0] << ", " << p.time_value[1] << ") "
- << "POL(" << pol_str << ")";
+ << "TIME[" << to_string(p.timesys) << "; " << p.time_value[0] << ", " << p.time_value[1] << "] "
+ << "POL[" << pol_str << "]";
return out;
}
--
GitLab