From 01ba8a157cf502a5993ed57eab980bfe80c069f4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pablo=20G=C3=B3mez?= Date: Wed, 8 Jul 2026 11:12:15 +0200 Subject: [PATCH 1/2] chore(release): v0.3.2 Hotfix release 0.3.2. Syncs the public mirror with the private v0.3.2 tag (excluding internal-only folders), covering: - FITS cutout WCS re-tangenting fix (parent-tile WCS fidelity) - UNIT/CONSVFLX FITS header keywords on individual cutouts - "Raw cutout (in Jy)" UI label and flux-conversion README docs - fitsbolt pinned to 0.3.0 --- CHANGELOG.md | 15 + README.md | 11 +- cutana/__init__.py | 2 +- cutana/cutout_extraction.py | 69 +++- cutana/cutout_process_utils.py | 9 + cutana/cutout_writer_fits.py | 232 +++++++----- cutana_ui/widgets/configuration_widget.py | 2 +- environment.yml | 2 +- pyproject.toml | 4 +- .../e2e/test_e2e_channel_combinations.py | 225 ++++++++---- tests/cutana/unit/test_cutout_extraction.py | 128 +++---- tests/cutana/unit/test_cutout_writer_fits.py | 330 +++++++++++++++++- 12 files changed, 788 insertions(+), 241 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index aa1d216..5d9c8a8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,21 @@ # Changelog +## [v0.3.2] – 2026-07-06 + +### Fixed +- **FITS cutout WCS re-tangenting bug**: cutout WCS now reproduces the parent tile mapping exactly — `CRVAL`, `CTYPE` and the CD/PC orientation are inherited from the parent unchanged and only `CRPIX` is shifted to the extraction origin — instead of re-tangenting the projection at each source. The previous approach left the cutout frame rotated by the meridian convergence between tile centre and source, producing a positional error that was ~0 at the cutout centre and grew toward the edges (~1″ at a few-arcmin FOV for sources far from the tile centre, worse near high \|Dec\| / tile corners). Only the FITS WCS header was affected; pixel data and Zarr/streaming outputs were not + +### Added +- **`UNIT` and `CONSVFLX` FITS header keywords** on individual cutout outputs: `UNIT` records the pixel unit (`OriginalUnit`, `Jy`, `approx Jy`, or `approx OriginalUnit`) depending on whether flux conversion and flux-conserved resizing were applied, and `CONSVFLX` records whether flux-conserved resizing was used + +### Changed +- **UI "Raw cutout" checkbox** relabelled to "Raw cutout (in Jy):" to make the output unit explicit +- **`fitsbolt` pinned to `==0.3.0`** (from `==0.2.0`) + +### Documentation +- **README**: added a Flux Conversion section (pixels converted to Jansky by default via the `MAGZERO` keyword, disable with `config.apply_flux_conversion = False`), documented the `UNIT` / `CONSVFLX` header keywords, and clarified the "Raw cutouts" terminology + ## [v0.3.1] – 2026-06-04 ### Fixed diff --git a/README.md b/README.md index 452b1da..099447b 100644 --- a/README.md +++ b/README.md @@ -17,7 +17,7 @@ ![Cutana Demo](assets/cutana_demo_2x.gif) > **Note:** Cutana is currently optimised for **Euclid Q1/IDR1 data**. Some defaults and assumptions are Euclid-specific: -> - Flux conversion expects the `MAGZERO` header keyword (configurable via `config.flux_conversion_keywords.AB_zeropoint`) +> - Flux conversion expects the `MAGZERO` header keyword (configurable via `config.flux_conversion_keywords.AB_zeropoint`) to convert to Jy > - Filter detection patterns are tuned for Euclid bands (VIS, NIR-Y, NIR-H, NIR-J) > - FITS structure assumes one file per channel/filter > @@ -100,8 +100,11 @@ TILE_102018666_12346,45.124,12.457,256,"['/path/to/tile_vis.fits','/path/to/tile **ZARR Format** (recommended): All cutouts stored in a efficient archives, ideal for large datasets and analysis workflows. Cutana uses the [Zarr format](https://zarr.readthedocs.io/en/stable/) for high-performance storage and the [images_to_zarr](https://github.com/gomezzz/images_to_zarr/) library for conversion. (See the Output section below for sample code to access) -**FITS Format**: Individual FITS files per source, best for compatibility with existing astronomical software. Mandatory format for `do_only_cutout_extraction`, -which skips all processing aside from the flux converison, which can be disabled. +**FITS Format**: Individual FITS files per source, best for compatibility with existing astronomical software. Mandatory format for `do_only_cutout_extraction` (Raw cutouts), +which skips all processing aside from the flux conversion, which can be disabled. + +### Flux Conversion +To align cutout units, pixels are converted to Janskys by default using the `MAGZERO` header keyword (configurable via `config.flux_conversion_keywords.AB_zeropoint`). Disable this flux conversion via `config.apply_flux_conversion = False`. ## WCS (World Coordinate System) Handling @@ -112,6 +115,8 @@ FITS cutouts **preserve full WCS information** with accurate astrometric calibra - **Reference coordinate centering**: WCS reference pixel (`CRPIX`) is set to the cutout center, with reference coordinates (`CRVAL`) pointing to the source position - **Format compatibility**: Supports CD matrix, CDELT, and PC+CDELT WCS formats from original FITS files - **Sky area preservation**: Total sky coverage remains constant while pixel scale adjusts for resize operations +- The header denotes the image unit (`UNIT`), `OriginalUnit`,`Jy`,`approx OriginalUnit`, `approx Jy`. +- `CONSVFLX` denotes if flux conserved resizing was applied. If not units are "approx" ### Zarr Output **Important**: Zarr archives **do not contain WCS information**. The WCS data is not recorded in the image metadata stored within the Zarr files. The central point and image size (in pixels or arcseconds, depending on what is provided) is recorded. diff --git a/cutana/__init__.py b/cutana/__init__.py index eb25053..f722023 100644 --- a/cutana/__init__.py +++ b/cutana/__init__.py @@ -27,7 +27,7 @@ # Application entry points (Orchestrator, UI) will enable logging when needed logger.disable("cutana") -__version__ = "0.3.1" +__version__ = "0.3.2" __author__ = "ESA Datalabs" # Import main classes for easy access diff --git a/cutana/cutout_extraction.py b/cutana/cutout_extraction.py index 1bec113..7eea1f7 100644 --- a/cutana/cutout_extraction.py +++ b/cutana/cutout_extraction.py @@ -106,11 +106,17 @@ def extract_cutouts_vectorized_from_extension( padding_factor: Factor to scale the extraction area (1.0 = no padding) Returns: - Tuple of (cutout_list, success_mask, pixel_offset_x, pixel_offset_y) where: + Tuple of (cutout_list, success_mask, pixel_offset_x, pixel_offset_y, + origin_x, origin_y) where: - cutout_list: List of cutout arrays (or None for failures) - success_mask: Boolean array indicating successful extractions - pixel_offset_x: Array of sub-pixel X offsets (positive = target toward right) - pixel_offset_y: Array of sub-pixel Y offsets (positive = target toward top) + - origin_x: Array of 0-based parent-pixel X origins of cutout pixel 0 (the + clip- and centre-pad-corrected window start). Pre-resize cutout pixel ``p`` + maps to parent pixel ``origin + p``; the FITS writer uses this to build the + cutout WCS without recomputing the (already-vectorised) extraction geometry. + - origin_y: As ``origin_x`` for the Y axis. """ n_sources = len(ra_array) logger.debug(f"Starting vectorized cutout extraction for {n_sources} sources") @@ -125,8 +131,10 @@ def extract_cutouts_vectorized_from_extension( return ( [None] * n_sources, np.zeros(n_sources, dtype=bool), - np.zeros(n_sources, dtype=np.float64), - np.zeros(n_sources, dtype=np.float64), + np.zeros(n_sources, dtype=np.float64), # sub-pixel offset x + np.zeros(n_sources, dtype=np.float64), # sub-pixel offset y + np.zeros(n_sources, dtype=np.int32), # integer origin x + np.zeros(n_sources, dtype=np.int32), # integer origin y ) img_height, img_width = image_data.shape @@ -147,8 +155,10 @@ def extract_cutouts_vectorized_from_extension( return ( [None] * n_sources, np.zeros(n_sources, dtype=bool), - np.zeros(n_sources, dtype=np.float64), - np.zeros(n_sources, dtype=np.float64), + np.zeros(n_sources, dtype=np.float64), # sub-pixel offset x + np.zeros(n_sources, dtype=np.float64), # sub-pixel offset y + np.zeros(n_sources, dtype=np.int32), # integer origin x + np.zeros(n_sources, dtype=np.int32), # integer origin y ) # Step 3: Vectorized bound computation @@ -191,6 +201,13 @@ def extract_cutouts_vectorized_from_extension( y_mins_clipped = np.maximum(0, y_mins) y_maxs_clipped = np.minimum(img_height, y_maxs) + # Parent-pixel origin of cutout pixel 0 (0-based, integer). For an on-tile window + # this is the clipped window start; edge-clipped windows are centre-padded below, + # which shifts the origin left/down by the integer pad offset. Threaded out so the + # FITS writer can build the cutout WCS without recomputing this geometry. + origin_x_array = x_mins_clipped.astype(np.int32) + origin_y_array = y_mins_clipped.astype(np.int32) + # Check for valid regions (vectorized) valid_mask = (x_maxs_clipped > x_mins_clipped) & (y_maxs_clipped > y_mins_clipped) @@ -263,6 +280,11 @@ def extract_cutouts_vectorized_from_extension( pixel_offset_x[i] -= pad_x_start pixel_offset_y[i] -= pad_y_start + # The data now starts pad pixels into the window, so cutout pixel 0 + # maps to a parent pixel pad_start before the clipped window start. + origin_x_array[i] -= pad_x_start + origin_y_array[i] -= pad_y_start + raw_cutout = padded_extraction # apply flux conversion here @@ -292,7 +314,7 @@ def extract_cutouts_vectorized_from_extension( successful_count = np.sum(success_mask) logger.debug(f"Vectorized extraction completed: {successful_count}/{n_sources} successful") - return cutouts, success_mask, pixel_offset_x, pixel_offset_y + return cutouts, success_mask, pixel_offset_x, pixel_offset_y, origin_x_array, origin_y_array def extract_cutouts_batch_vectorized( @@ -382,17 +404,22 @@ def extract_cutouts_batch_vectorized( logger.debug(f"Processing extension {ext_name} for {n_sources} sources") # Extract cutouts for all sources in this extension using vectorized method - cutout_list, success_mask, offset_x_array, offset_y_array = ( - extract_cutouts_vectorized_from_extension( - hdul[ext_name], - wcs_dict[ext_name], - ra_array, - dec_array, - size_pixels_array, - source_ids, - padding_factor, - config, - ) + ( + cutout_list, + success_mask, + offset_x_array, + offset_y_array, + origin_x_array, + origin_y_array, + ) = extract_cutouts_vectorized_from_extension( + hdul[ext_name], + wcs_dict[ext_name], + ra_array, + dec_array, + size_pixels_array, + source_ids, + padding_factor, + config, ) # Organize results by source ID @@ -401,10 +428,16 @@ def extract_cutouts_batch_vectorized( if source_id not in combined_cutouts: combined_cutouts[source_id] = {} combined_wcs[source_id] = {} - # Store pixel offsets (same for all extensions since coords are the same) + # Store per-source geometry (same for all extensions since the + # coordinates and window are identical): the sub-pixel offsets and + # the integer extraction origin/size the FITS writer needs to build + # the cutout WCS without recomputing world_to_pixel and the bounds. combined_offsets[source_id] = { "x": float(offset_x_array[i]), "y": float(offset_y_array[i]), + "origin_x": int(origin_x_array[i]), + "origin_y": int(origin_y_array[i]), + "extraction_size": int(size_pixels_array[i] * padding_factor), } combined_cutouts[source_id][ext_name] = cutout diff --git a/cutana/cutout_process_utils.py b/cutana/cutout_process_utils.py index f78a2d2..b6979cb 100644 --- a/cutana/cutout_process_utils.py +++ b/cutana/cutout_process_utils.py @@ -493,6 +493,12 @@ def _process_sources_batch_vectorized_with_fits_set( source_offsets = all_source_offsets.get(source_id, {"x": 0.0, "y": 0.0}) extraction_offset_x = source_offsets.get("x", 0.0) extraction_offset_y = source_offsets.get("y", 0.0) + # Integer extraction origin/size (parent-tile pixels, pre-resize), computed + # vectorised at extraction time. Passed through so the FITS writer can build + # the cutout WCS without recomputing world_to_pixel and the window bounds. + extraction_origin_x = source_offsets.get("origin_x") + extraction_origin_y = source_offsets.get("origin_y") + extraction_size = source_offsets.get("extraction_size") # When resizing is applied, scale offsets and pixel scale by the same # resize_factor so both stay consistent with the final output coords. @@ -541,6 +547,9 @@ def _process_sources_batch_vectorized_with_fits_set( "processing_timestamp": batch_timestamp, "rescaled_offset_x": rescaled_offset_x, "rescaled_offset_y": rescaled_offset_y, + "extraction_origin_x": extraction_origin_x, + "extraction_origin_y": extraction_origin_y, + "extraction_size": extraction_size, } ) wcs_list.append(all_source_wcs.get(source_id, {})) diff --git a/cutana/cutout_writer_fits.py b/cutana/cutout_writer_fits.py index 6037a79..22c28b4 100644 --- a/cutana/cutout_writer_fits.py +++ b/cutana/cutout_writer_fits.py @@ -114,36 +114,87 @@ def generate_fits_filename( return f"{source_id}_cutout.fits" +def _rescale_wcs_pixel_scale(header: fits.Header, pixel_scale_matrix: Any, resize: float) -> None: + """Rescale a WCS header's pixel scale (CD or PC+CDELT) in place for a resize. + + A resize by ``resize = final_size / extraction_size`` shrinks the per-pixel sky + step by ``1 / resize`` (more output pixels -> smaller step). The projection + orientation (CRVAL, CTYPE, and the CD/PC direction) is untouched. + + Args: + header: WCS header to modify in place. + pixel_scale_matrix: Parent ``pixel_scale_matrix`` (or None if unavailable). + resize: Ratio final_size / extraction_size. + """ + if "CD1_1" in header and "CD2_2" in header: + # CD already folds scale into the matrix; divide every present element. + for key in ("CD1_1", "CD1_2", "CD2_1", "CD2_2"): + if key in header: + header[key] = header[key] / resize + elif "CDELT1" in header and "CDELT2" in header: + if pixel_scale_matrix is not None and "PC1_1" in header and "PC2_2" in header: + # PC carries the orientation; put the rescaled per-pixel step into CDELT. + header["CDELT1"] = (pixel_scale_matrix[0, 0] / resize) / header.get("PC1_1", 1.0) + header["CDELT2"] = (pixel_scale_matrix[1, 1] / resize) / header.get("PC2_2", 1.0) + else: + header["CDELT1"] = header["CDELT1"] / resize + header["CDELT2"] = header["CDELT2"] / resize + + def create_wcs_header( cutout_shape: tuple, original_wcs: Optional[WCS] = None, - ra_center: Optional[float] = None, - dec_center: Optional[float] = None, + ra_source: Optional[float] = None, + dec_source: Optional[float] = None, pixel_scale: Optional[float] = None, resize_factor: Optional[float] = None, rescaled_offset_x: Optional[float] = None, rescaled_offset_y: Optional[float] = None, + extraction_origin_x: Optional[int] = None, + extraction_origin_y: Optional[int] = None, + extraction_size: Optional[int] = None, ) -> fits.Header: """ - Create WCS header for cutout. + Create a WCS header for a cutout. + + When the parent-tile WCS is available, the cutout WCS reproduces the parent + sky mapping exactly (the ``astropy.nddata.Cutout2D`` construction): CRVAL, + CTYPE and the CD/PC orientation are inherited unchanged, and only CRPIX is + shifted to the cutout's integer extraction origin (CD rescaled when the + cutout was resized). + + It deliberately does NOT re-tangent the projection at the source position. + Setting CRVAL to the source RA/Dec and CRPIX to the geometric centre while + keeping the tile's CD matrix rotates the cutout frame by the meridian + convergence between the tile centre and the source, producing a WCS error + that GROWS with distance from the cutout centre (order ~1" at a few-arcmin + FOV for sources far from the tile centre). + + The extraction origin and size are threaded in from ``cutout_extraction`` (via + metadata) rather than recomputed here, so the (already-vectorised) world_to_pixel + and window/clip/pad geometry is not duplicated. Args: - cutout_shape: Shape of the cutout (height, width) - original_wcs: Original WCS from parent image - ra_center: RA of cutout center in degrees - dec_center: Dec of cutout center in degrees - pixel_scale: Pixel scale in arcsec/pixel - resize_factor: Factor by which the cutout was resized (new_size/original_size) - Used ONLY for adjusting pixel scale in WCS, NOT for offset scaling. - rescaled_offset_x: Sub-pixel X offset in FINAL image coordinates (positive = target toward right). - This offset is ALREADY scaled by resize_factor and should be used as-is. - rescaled_offset_y: Sub-pixel Y offset in FINAL image coordinates (positive = target toward top). - This offset is ALREADY scaled by resize_factor and should be used as-is. + cutout_shape: Shape of the final cutout (height, width). + original_wcs: WCS of the parent tile. + ra_source: Source RA in degrees. Used as CRVAL only in the no-parent-WCS fallback. + dec_source: Source Dec in degrees (see ``ra_source``). + pixel_scale: Pixel scale in arcsec/pixel (fallback branch only). + resize_factor: final_size / requested_size (fallback branch only). + rescaled_offset_x: Sub-pixel X offset in final coords (fallback branch only). + rescaled_offset_y: Sub-pixel Y offset in final coords (fallback branch only). + extraction_origin_x: 0-based integer parent-pixel X origin of cutout pixel 0, + computed at extraction time (clip/centre-pad corrected). Required for the + parent-WCS branch. + extraction_origin_y: As ``extraction_origin_x`` for the Y axis. + extraction_size: Pre-resize extraction window size in parent pixels + (``int(requested * padding_factor)``). Required for the parent-WCS branch; + the resize ratio is ``final_size / extraction_size``. Returns: - FITS header with WCS information + FITS header with WCS information. """ - # Default offsets to 0 if not provided + # Default offsets to 0 if not provided (fallback branch only). if rescaled_offset_x is None: rescaled_offset_x = 0.0 if rescaled_offset_y is None: @@ -155,27 +206,53 @@ def create_wcs_header( cached_header, cached_pixel_scale_matrix = _get_cached_wcs_info(original_wcs) header = cached_header.copy() - # Update reference pixel to center of cutout, adjusted by rescaled offset - # CRPIX follows FITS convention: 1-based indexing where pixel (1,1) is bottom-left - # For an N-pixel image, the geometric center is at (N/2 + 0.5) in FITS 1-based coords - # The rescaled_offset is in 0-based pixel coordinates, so we add it after converting center to 1-based - height, width = cutout_shape - fits_center_x = width / 2.0 + 0.5 # Convert 0-based center to FITS 1-based - fits_center_y = height / 2.0 + 0.5 - header["CRPIX1"] = fits_center_x + rescaled_offset_x - header["CRPIX2"] = fits_center_y + rescaled_offset_y + final_size = cutout_shape[0] # cutouts are square + + if ( + extraction_origin_x is None + or extraction_origin_y is None + or extraction_size is None + ): + # These are threaded from cutout_extraction via metadata; their absence + # is a broken invariant, not a recoverable per-source condition. + raise ValueError( + "extraction_origin_x/y and extraction_size are required to build a " + "cutout WCS from a parent tile WCS (threaded from cutout_extraction)." + ) + + if getattr(original_wcs, "sip", None) is not None: + # SIP distortion is referenced to CRPIX and is not propagated/rescaled + # here, so a distorted tile would yield an inaccurate cutout WCS. Euclid + # MER mosaics are undistorted TAN, so this does not affect them; warn + # loudly rather than emit a silently-wrong header for distorted tiles. + # Full SIP propagation is tracked in issue #238. + logger.warning( + "Parent tile WCS carries SIP distortion, which the cutout WCS does " + "not propagate; the cutout WCS may be inaccurate for distorted tiles " + "(full SIP support tracked in issue #238)." + ) + + resize = final_size / extraction_size + crpix1_parent = float(cached_header["CRPIX1"]) + crpix2_parent = float(cached_header["CRPIX2"]) + # Pre-resize: CRPIX_cut = CRPIX_parent - origin. A resize by ``resize`` maps + # this to 0.5 + resize*(CRPIX - origin - 0.5) (cv2.resize half-pixel centres). + header["CRPIX1"] = 0.5 + resize * (crpix1_parent - extraction_origin_x - 0.5) + header["CRPIX2"] = 0.5 + resize * (crpix2_parent - extraction_origin_y - 0.5) logger.debug( - f"WCS CRPIX: FITS_center=({fits_center_x:.2f}, {fits_center_y:.2f}) + " - f"offset=({rescaled_offset_x:.4f}, {rescaled_offset_y:.4f}) = " + f"WCS CRPIX (origin-anchored): ext={extraction_size} resize={resize:.4f} " + f"origin=({extraction_origin_x}, {extraction_origin_y}) -> " f"({header['CRPIX1']:.4f}, {header['CRPIX2']:.4f})" ) - # Update reference coordinates if provided - if ra_center is not None and dec_center is not None: - header["CRVAL1"] = ra_center - header["CRVAL2"] = dec_center + # CRVAL, CTYPE and the CD/PC orientation are inherited from the parent + # unchanged (do NOT overwrite CRVAL). Only rescale the pixel scale for a resize. + if resize != 1.0: + _rescale_wcs_pixel_scale(header, cached_pixel_scale_matrix, resize) + + return header - elif ra_center is not None and dec_center is not None: + elif ra_source is not None and dec_source is not None: # Create minimal WCS header header = fits.Header() height, width = cutout_shape @@ -188,8 +265,8 @@ def create_wcs_header( fits_center_y = height / 2.0 + 0.5 header["CRPIX1"] = fits_center_x + rescaled_offset_x header["CRPIX2"] = fits_center_y + rescaled_offset_y - header["CRVAL1"] = ra_center - header["CRVAL2"] = dec_center + header["CRVAL1"] = ra_source + header["CRVAL2"] = dec_source logger.debug( f"Minimal WCS CRPIX: FITS_center=({fits_center_x:.2f}, {fits_center_y:.2f}) + " f"offset=({rescaled_offset_x:.4f}, {rescaled_offset_y:.4f}) = " @@ -221,52 +298,6 @@ def create_wcs_header( # No WCS info available return fits.Header() - # Apply resize factor to pixel scale when we have original_wcs - if original_wcs is not None and resize_factor is not None and resize_factor != 1.0: - # Scale pixel scale by the resize factor - # If image was made smaller (resize_factor < 1), pixels represent larger sky area - # If image was made larger (resize_factor > 1), pixels represent smaller sky area - - # Use cached pixel scale matrix (computed once per WCS) - if cached_pixel_scale_matrix is not None: - original_pixel_scale_x = cached_pixel_scale_matrix[0, 0] - original_pixel_scale_y = cached_pixel_scale_matrix[1, 1] - - # Apply resize factor to get new pixel scale - new_pixel_scale_x = original_pixel_scale_x / resize_factor - new_pixel_scale_y = original_pixel_scale_y / resize_factor - - # Handle CD matrix (preferred modern format) - if "CD1_1" in header and "CD2_2" in header: - header["CD1_1"] = header["CD1_1"] / resize_factor - header["CD2_2"] = header["CD2_2"] / resize_factor - if "CD1_2" in header: - header["CD1_2"] = header["CD1_2"] / resize_factor - if "CD2_1" in header: - header["CD2_1"] = header["CD2_1"] / resize_factor - - # Handle CDELT format or PC+CDELT format - elif "CDELT1" in header and "CDELT2" in header: - # For PC+CDELT format, set CDELT to achieve desired pixel scale - if "PC1_1" in header and "PC2_2" in header: - pc1_1 = header.get("PC1_1", 1.0) - pc2_2 = header.get("PC2_2", 1.0) - header["CDELT1"] = new_pixel_scale_x / pc1_1 - header["CDELT2"] = new_pixel_scale_y / pc2_2 - else: - header["CDELT1"] = new_pixel_scale_x - header["CDELT2"] = new_pixel_scale_y - else: - # Fallback: simple scaling of existing header values - if "CD1_1" in header and "CD2_2" in header: - header["CD1_1"] = header["CD1_1"] / resize_factor - header["CD2_2"] = header["CD2_2"] / resize_factor - elif "CDELT1" in header and "CDELT2" in header: - header["CDELT1"] = header["CDELT1"] / resize_factor - header["CDELT2"] = header["CDELT2"] / resize_factor - - return header - except Exception as e: logger.error(f"Failed to create WCS header: {e}") # Return minimal header @@ -327,6 +358,11 @@ def write_single_fits_cutout( "DTYPE": metadata.get("data_type", "float32"), "PIXSCALE": metadata.get("pixel_scale_arcsec_per_pixel"), "TILE": metadata.get("tile"), + "UNIT": (cutout_data.get("unit", "---"), "Data unit, 'approx' = non-flux conserved"), + "CONSVFLX": ( + cutout_data.get("conserved_flux", False), + "F: Flux not conserved in resizing, T: conserved", + ), } primary_hdu.header.update(primary_header_updates) @@ -366,6 +402,13 @@ def write_single_fits_cutout( f"Retrieved rescaled offsets from metadata: ({rescaled_offset_x:.4f}, {rescaled_offset_y:.4f})" ) + # Integer extraction origin/size threaded from cutout_extraction; the + # parent-WCS branch uses these to anchor CRPIX without recomputing the + # window geometry. + extraction_origin_x = metadata.get("extraction_origin_x") + extraction_origin_y = metadata.get("extraction_origin_y") + extraction_size = metadata.get("extraction_size") + if channel in wcs_info: logger.debug( f"Creating WCS header for channel {channel} using original WCS" @@ -373,11 +416,14 @@ def write_single_fits_cutout( wcs_header = create_wcs_header( cutout.shape, original_wcs=wcs_info[channel], - ra_center=metadata.get("ra"), - dec_center=metadata.get("dec"), + ra_source=metadata.get("ra"), + dec_source=metadata.get("dec"), resize_factor=resize_factor, rescaled_offset_x=rescaled_offset_x, rescaled_offset_y=rescaled_offset_y, + extraction_origin_x=extraction_origin_x, + extraction_origin_y=extraction_origin_y, + extraction_size=extraction_size, ) else: # Fallback: create minimal WCS using source coordinates @@ -387,11 +433,14 @@ def write_single_fits_cutout( wcs_header = create_wcs_header( cutout.shape, original_wcs=None, - ra_center=metadata.get("ra"), - dec_center=metadata.get("dec"), + ra_source=metadata.get("ra"), + dec_source=metadata.get("dec"), resize_factor=resize_factor, rescaled_offset_x=rescaled_offset_x, rescaled_offset_y=rescaled_offset_y, + extraction_origin_x=extraction_origin_x, + extraction_origin_y=extraction_origin_y, + extraction_size=extraction_size, ) if wcs_header: @@ -522,8 +571,23 @@ def write_fits_batch( "metadata": metadata, "processed_cutouts": processed_cutouts, "wcs_info": source_wcs_info, # Use properly mapped WCS info + "unit": "OriginalUnit", + "conserved_flux": False, } + if config.flux_conserved_resizing or config.do_only_cutout_extraction: + cutout_data["conserved_flux"] = True + if config.apply_flux_conversion: + cutout_data["unit"] = ( + "Jy" # Example unit after flux conversion, adjust as needed + ) + + else: + if config.apply_flux_conversion: + cutout_data["unit"] = "approx Jy" + else: + cutout_data["unit"] = "approx OriginalUnit" + # Determine output directory for this source if create_subdirs: # Create subdirectory based on first few characters of source ID diff --git a/cutana_ui/widgets/configuration_widget.py b/cutana_ui/widgets/configuration_widget.py index aee0611..3771d98 100644 --- a/cutana_ui/widgets/configuration_widget.py +++ b/cutana_ui/widgets/configuration_widget.py @@ -149,7 +149,7 @@ def __init__( # Raw cutout only checkbox - disables all processing when checked self.do_only_cutout_label = widgets.HTML( - value=f'
Raw cutout:
', + value=f'
Raw cutout [Jy]:
', layout=widgets.Layout(height="28px", width="100%"), ) self.do_only_cutout_checkbox = widgets.Checkbox( diff --git a/environment.yml b/environment.yml index 609af4c..7205dfc 100644 --- a/environment.yml +++ b/environment.yml @@ -37,7 +37,7 @@ dependencies: - vulture>=2.10 # Dead code detection - pip: - drizzle>=2.0.1 # Image resampling with drizzle algorithm - - fitsbolt==0.2.0 # FITS file loading and normalisation + - fitsbolt==0.3.0 # FITS file loading and normalisation - images-to-zarr>=0.3.5 # Image to Zarr conversion - ipyfilechooser>=0.6 # File/folder selection widget - memory-profiler>=0.60 # Memory usage profiling diff --git a/pyproject.toml b/pyproject.toml index 951e9a5..408d60d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "cutana" -version = "0.3.1" +version = "0.3.2" description = "High-performance Python tool for creating astronomical image cutouts" readme = "README.md" authors = [{ name = "European Space Agency", email = "pablo.gomez@esa.int" }] @@ -20,7 +20,7 @@ dependencies = [ "astropy>=5.0", "dotmap>=1.3", "drizzle>=2.0.1", - "fitsbolt==0.2.0", + "fitsbolt==0.3.0", "fsspec>=2023.1.0", "images-to-zarr>=0.3.5", "ipyfilechooser>=0.6", diff --git a/tests/cutana/e2e/test_e2e_channel_combinations.py b/tests/cutana/e2e/test_e2e_channel_combinations.py index 9b8c5eb..e7d67ef 100644 --- a/tests/cutana/e2e/test_e2e_channel_combinations.py +++ b/tests/cutana/e2e/test_e2e_channel_combinations.py @@ -1556,9 +1556,11 @@ def test_fits_output_wcs_preservation(self, temp_dir, mock_fits_files, mode, flu This test validates that: 1. Output FITS files have valid WCS headers - 2. The WCS reference coordinates (CRVAL) match the source RA/Dec + 2. CRVAL/CTYPE are inherited from the parent tile (projection not re-tangented; + in this fixture the tile is centred on the source, so CRVAL == source RA/Dec) 3. The pixel scale is correctly adjusted for resizing (original -> target resolution) - 4. The WCS can correctly convert pixel coordinates to sky coordinates + 4. The output WCS reproduces the parent-tile mapping: the source RA/Dec lands at + the pixel where the extracted data actually sits, not the geometric centre Test modes: - standard: Normal processing with resize (default behavior) @@ -1760,79 +1762,50 @@ def test_fits_output_wcs_preservation(self, temp_dir, mock_fits_files, mode, flu f"CRVAL2 (Dec) mismatch: expected {test_dec}, got {crval_dec}" ) - # --- Independently calculate the expected pixel offset --- - + # --- Verify the output WCS reproduces the parent-tile mapping --- + # The cutout WCS must locate the source RA/Dec where the extraction + # (integer window + cv2 half-pixel-centre resize) actually places the + # data, NOT at the geometric centre. Ground truth is derived from the + # parent tile WCS and the extraction geometry, independent of the WCS + # construction under test. orig_wcs = mock_wcs_data skycoord = SkyCoord(ra=test_ra * u.deg, dec=test_dec * u.deg, frame="icrs") pixel_x, pixel_y = orig_wcs.world_to_pixel(skycoord) - extraction_size = 10 + extraction_size = original_cutout_size half_left = extraction_size // 2 - # Use np.floor to match .astype(int) behavior in main code - x_min = (np.asarray(pixel_x - half_left)).astype(int) - y_min = (np.asarray(pixel_y - half_left)).astype(int) - cutout_center_x = x_min + extraction_size / 2.0 - cutout_center_y = y_min + extraction_size / 2.0 - pixel_offset_x = pixel_x - cutout_center_x - pixel_offset_y = pixel_y - cutout_center_y - + x_min = int(np.asarray(pixel_x - half_left)) + y_min = int(np.asarray(pixel_y - half_left)) resize_factor = expected_resize_factor - if resize_factor != 1.0: - pixel_offset_x *= resize_factor - pixel_offset_y *= resize_factor - # Calculate expected CRPIX using FITS 1-based indexing - # For an N-pixel image, center is at (N/2 + 0.5) in FITS 1-based coordinates - fits_center_x = cutout_shape[1] / 2.0 + 0.5 - fits_center_y = cutout_shape[0] / 2.0 + 0.5 - expected_crpix1 = fits_center_x + pixel_offset_x - expected_crpix2 = fits_center_y + pixel_offset_y + # Expected 0-based source position in the (resized) cutout, using the + # cv2.resize half-pixel-centre convention: out = (in + 0.5) * r - 0.5. + exp_x0 = (float(pixel_x) - x_min + 0.5) * resize_factor - 0.5 + exp_y0 = (float(pixel_y) - y_min + 0.5) * resize_factor - 0.5 - crpix1 = wcs.wcs.crpix[0] - crpix2 = wcs.wcs.crpix[1] - - logger.info(f" WCS CRPIX (FITS 1-based): ({crpix1}, {crpix2})") + out_x0, out_y0 = wcs.world_to_pixel(skycoord) logger.info( - f" Expected CRPIX (FITS 1-based): ({expected_crpix1}, {expected_crpix2})" + f" Source output pixel (0-based): got=({out_x0}, {out_y0}) " + f"expected=({exp_x0}, {exp_y0})" ) - # Get the pixel position for the original RA/Dec using the output WCS - # world_to_pixel returns 0-based pixel coordinates - pixel_from_wcs_0based = wcs.world_to_pixel(skycoord) - # Convert to FITS 1-based for comparison - pixel_from_wcs_1based_x = pixel_from_wcs_0based[0] + 1 - pixel_from_wcs_1based_y = pixel_from_wcs_0based[1] + 1 - logger.info(f" Output WCS pixel for RA/Dec (0-based): {pixel_from_wcs_0based}") - logger.info( - f" Output WCS pixel for RA/Dec (FITS 1-based): ({pixel_from_wcs_1based_x}, {pixel_from_wcs_1based_y})" + # The fix reproduces the parent mapping exactly, so the output + # pixel must match the extraction-derived expectation to well within + # a hundredth of a pixel (this is not a discretisation tolerance). + pixel_tolerance = 0.05 + assert abs(float(out_x0) - exp_x0) < pixel_tolerance, ( + f"Source X pixel mismatch: expected {exp_x0}, got {out_x0}" ) - - # All comparisons use FITS 1-based indexing - # Tolerance: 0.05 * original pixel size (in pixels, not degrees) - pixel_tolerance = 0.05 * original_cutout_size - - # Compare CRPIX to expected (both FITS 1-based) - assert abs(crpix1 - expected_crpix1) < pixel_tolerance, ( - f"CRPIX1 mismatch: expected {expected_crpix1}, got {crpix1}, tol={pixel_tolerance}" - ) - assert abs(crpix2 - expected_crpix2) < pixel_tolerance, ( - f"CRPIX2 mismatch: expected {expected_crpix2}, got {crpix2}, tol={pixel_tolerance}" - ) - - # Compare output WCS pixel (converted to FITS 1-based) to expected CRPIX (FITS 1-based) - assert abs(pixel_from_wcs_1based_x - expected_crpix1) < pixel_tolerance, ( - f"Output WCS pixel X for RA/Dec mismatch: expected {expected_crpix1}, got {pixel_from_wcs_1based_x}, tol={pixel_tolerance}" - ) - assert abs(pixel_from_wcs_1based_y - expected_crpix2) < pixel_tolerance, ( - f"Output WCS pixel Y for RA/Dec mismatch: expected {expected_crpix2}, got {pixel_from_wcs_1based_y}, tol={pixel_tolerance}" + assert abs(float(out_y0) - exp_y0) < pixel_tolerance, ( + f"Source Y pixel mismatch: expected {exp_y0}, got {out_y0}" ) - # Compare output WCS pixel to CRPIX (both FITS 1-based) - assert abs(pixel_from_wcs_1based_x - crpix1) < pixel_tolerance, ( - f"Output WCS pixel X for RA/Dec mismatch with CRPIX1: {pixel_from_wcs_1based_x} vs {crpix1}, tol={pixel_tolerance}" - ) - assert abs(pixel_from_wcs_1based_y - crpix2) < pixel_tolerance, ( - f"Output WCS pixel Y for RA/Dec mismatch with CRPIX2: {pixel_from_wcs_1based_y} vs {crpix2}, tol={pixel_tolerance}" + # Cross-check: at the pixel where the source data actually sits, the + # output WCS must report the source's true sky position (sub-mas). + sky_at_source = wcs.pixel_to_world(exp_x0, exp_y0) + sep_arcsec = sky_at_source.separation(skycoord).to(u.arcsec).value + assert sep_arcsec < 1e-3, ( + f"Output WCS sky position at the source pixel is off by {sep_arcsec} arcsec" ) # Pixel scale check (unchanged) @@ -1865,6 +1838,136 @@ def test_fits_output_wcs_preservation(self, temp_dir, mock_fits_files, mode, flu f"WCS preservation test passed for mode={mode}, flux_conserved={flux_conserved}!" ) + def test_fits_output_wcs_offcenter_source_near_pole(self, temp_dir): + """E2E: a source far off the cutout centre in a large near-pole cutout. + + This is the scenario the re-tangenting bug hit hardest: a cutout whose centre + is 2 deg from the tile centre, with a source 30 arcmin from the cutout centre, + near the pole where meridian convergence is strong. A coarse 1 arcmin/px scale + keeps the mock tile tiny (400x400) for fast compute while still spanning ~6.7 deg. + + A single bright pixel marks the source's true sky position. After the full + pipeline, feeding that source's RA/Dec through the OUTPUT WCS must land on the + bright pixel (and the reverse: the bright pixel's sky position must match the + source). The old re-tangented WCS would miss it by arcsec-to-arcmin here. + """ + # --- coarse, large, near-pole tile --- + tile_npix = 400 + pixscale_deg = 60.0 / 3600.0 # 1 arcmin/px -> 400 px spans ~6.7 deg + tile_wcs = WCS(naxis=2) + tile_wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"] + tile_wcs.wcs.crpix = [tile_npix / 2.0, tile_npix / 2.0] + tile_wcs.wcs.crval = [30.0, -80.0] # near the south pole + tile_wcs.wcs.cd = [[-pixscale_deg, 0.0], [0.0, pixscale_deg]] + tile_wcs.wcs.cunit = ["deg", "deg"] + tile_wcs.pixel_shape = (tile_npix, tile_npix) + + # Cutout centre 2 deg from the tile centre (moved along +Dec, away from pole). + offset_px = int(round(2.0 / pixscale_deg)) # 120 px + target_px = tile_npix / 2.0 + target_py = tile_npix / 2.0 - offset_px + target_sky = tile_wcs.pixel_to_world(target_px, target_py) + target_ra, target_dec = float(target_sky.ra.deg), float(target_sky.dec.deg) + + # Bright 1-px source ~30 arcmin from the cutout centre, offset in BOTH axes + # (24' in x, 18' in y -> 30' radial) so the WCS is checked off both axes. + src_x_offset_px = int(round((24.0 / 60.0) / pixscale_deg)) # 24 px + src_y_offset_px = int(round((18.0 / 60.0) / pixscale_deg)) # 18 px + src_px = int(round(target_px)) + src_x_offset_px + src_py = int(round(target_py)) + src_y_offset_px + src_sky = tile_wcs.pixel_to_world(src_px, src_py) + src_coord = SkyCoord(ra=src_sky.ra.deg * u.deg, dec=src_sky.dec.deg * u.deg, frame="icrs") + + data = np.zeros((tile_npix, tile_npix), dtype=np.float32) + data[src_py, src_px] = 1000.0 + + # --- write tile to disk (for catalogue analysis) + mock the loader --- + filename = "EUC_MER_BGSUB-MOSAIC-VIS_TILE111111111-ABCDE_20001021T024200.000000Z_00.00.fits" + fits_path = Path(temp_dir) / filename + header = tile_wcs.to_header() + header["MAGZERO"] = 25.0 + header["EXTNAME"] = "PRIMARY" + header["BUNIT"] = "electron/s" + header["INSTRUME"] = "VIS" + fits.PrimaryHDU(data=data, header=header).writeto(fits_path, overwrite=True) + + # Cutout diameter large enough to contain the 30 px-offset source (radius 40 px). + diameter = 80 + catalogue_data = [ + { + "SourceID": "offcenter_src", + "RA": target_ra, + "Dec": target_dec, + "diameter_pixel": diameter, + "fits_file_paths": json.dumps([str(fits_path)]), + } + ] + catalogue_path = Path(temp_dir) / "offcenter_cat.csv" + pd.DataFrame(catalogue_data).to_csv(catalogue_path, index=False) + result_cat_analysis = analyse_source_catalogue(str(catalogue_path)) + + selected_fits = {"VIS": {"path": str(fits_path), "value": 100.0, "filename": filename}} + + config = get_default_config() + config.source_catalogue = str(catalogue_path) + config.output_dir = str(Path(temp_dir) / "offcenter_output") + config.data_type = "float32" + config.normalisation_method = "none" # keep the bright pixel intact + config.apply_flux_conversion = False + config.max_workers = 1 + config.output_format = "fits" + config.N_batch_cutout_process = 10 + config.do_only_cutout_extraction = True # raw 1:1 extraction, no resize/mixing + config.target_resolution = 16 # unused in cutout_only but must pass validation + config.flux_conserved_resizing = False + config.fits_extensions = ["PRIMARY"] + config.channel_weights, config.selected_extensions, config.available_extensions = ( + self.create_ui_shared_config_like_config(selected_fits, 1, result_cat_analysis, 1, 1) + ) + Path(config.output_dir).mkdir(parents=True, exist_ok=True) + + orchestrator = Orchestrator(config) + mock_fits_data = None + with patch("cutana.fits_dataset.load_fits_sets") as mock_load_fits: + hdul = fits.HDUList([fits.PrimaryHDU(data=data, header=header)]) + mock_fits_data = {str(fits_path): (hdul, {"PRIMARY": tile_wcs})} + mock_load_fits.return_value = mock_fits_data + + cat_parquet = Path(config.output_dir) / "cat.parquet" + pd.DataFrame(catalogue_data).to_parquet(cat_parquet, index=False) + config.source_catalogue = str(cat_parquet) + try: + result = orchestrator.start_processing(str(cat_parquet)) + assert result["status"] == "completed", f"Processing failed: {result}" + finally: + try: + orchestrator.stop_processing() + except Exception: + pass + self._close_mock_fits_data(mock_fits_data) + + # --- verify the source lands correctly under the OUTPUT WCS --- + output_files = list(Path(config.output_dir).glob("*.fits")) + assert len(output_files) == 1, f"expected one cutout, got {output_files}" + with fits.open(output_files[0]) as hdul_out: + data_ext = [h for h in hdul_out if getattr(h, "data", None) is not None][-1] + out_wcs = WCS(data_ext.header) + cut = data_ext.data + assert cut.shape == (diameter, diameter) + + # Bright pixel = the source in the extracted cutout. + cj, ci = np.unravel_index(int(np.argmax(cut)), cut.shape) # (row=y, col=x) + + # Forward: the source's true RA/Dec, via the output WCS, must land on it. + wx, wy = out_wcs.world_to_pixel(src_coord) + assert abs(float(wx) - ci) < 0.05, f"source X: WCS {wx} vs data {ci}" + assert abs(float(wy) - cj) < 0.05, f"source Y: WCS {wy} vs data {cj}" + + # Reverse: the bright pixel's sky position must match the source (sub-arcsec). + sky = out_wcs.pixel_to_world(int(ci), int(cj)) + sep = sky.separation(src_coord).to(u.arcsec).value + assert sep < 0.05, f"source WCS off by {sep} arcsec at ({ci}, {cj})" + def test_fits_output_wcs_combined_channels(self, temp_dir, mock_fits_files): """Test that WCS is correctly preserved when combining VIS and NIR-H into one output channel. diff --git a/tests/cutana/unit/test_cutout_extraction.py b/tests/cutana/unit/test_cutout_extraction.py index 1180dcf..1d98ac4 100644 --- a/tests/cutana/unit/test_cutout_extraction.py +++ b/tests/cutana/unit/test_cutout_extraction.py @@ -71,15 +71,17 @@ def test_odd_size_extraction_no_padding(self, mock_hdu_ones, mock_wcs): ra, dec = 150.0, 2.0 # These should map to approximately center for size in odd_sizes: - cutouts, success_mask, offset_x, offset_y = extract_cutouts_vectorized_from_extension( - mock_hdu_ones, - mock_wcs, - np.array([ra]), - np.array([dec]), - np.array([size]), - source_ids=["test_source"], - padding_factor=1.0, - config=None, + cutouts, success_mask, offset_x, offset_y, *_ = ( + extract_cutouts_vectorized_from_extension( + mock_hdu_ones, + mock_wcs, + np.array([ra]), + np.array([dec]), + np.array([size]), + source_ids=["test_source"], + padding_factor=1.0, + config=None, + ) ) assert success_mask[0], f"Extraction failed for size {size}" @@ -104,15 +106,17 @@ def test_even_size_extraction_no_padding(self, mock_hdu_ones, mock_wcs): ra, dec = 150.0, 2.0 for size in even_sizes: - cutouts, success_mask, offset_x, offset_y = extract_cutouts_vectorized_from_extension( - mock_hdu_ones, - mock_wcs, - np.array([ra]), - np.array([dec]), - np.array([size]), - source_ids=["test_source"], - padding_factor=1.0, - config=None, + cutouts, success_mask, offset_x, offset_y, *_ = ( + extract_cutouts_vectorized_from_extension( + mock_hdu_ones, + mock_wcs, + np.array([ra]), + np.array([dec]), + np.array([size]), + source_ids=["test_source"], + padding_factor=1.0, + config=None, + ) ) assert success_mask[0], f"Extraction failed for size {size}" @@ -138,15 +142,17 @@ def test_edge_extraction_with_padding(self, mock_hdu_ones, mock_wcs): mock_world_to_pixel.return_value = (5.0, 5.0) # Near top-left corner size = 20 - cutouts, success_mask, offset_x, offset_y = extract_cutouts_vectorized_from_extension( - mock_hdu_ones, - mock_wcs, - np.array([150.0]), - np.array([2.0]), - np.array([size]), - source_ids=["edge_source"], - padding_factor=1.0, - config=None, + cutouts, success_mask, offset_x, offset_y, *_ = ( + extract_cutouts_vectorized_from_extension( + mock_hdu_ones, + mock_wcs, + np.array([150.0]), + np.array([2.0]), + np.array([size]), + source_ids=["edge_source"], + padding_factor=1.0, + config=None, + ) ) assert success_mask[0], "Extraction failed for edge source" @@ -172,15 +178,17 @@ def test_flux_conversion_applied(self, mock_hdu_ones, mock_wcs): # Make flux conversion multiply by 2 for testing mock_flux_conv.return_value = np.ones((10, 10)) * 2.0 - cutouts, success_mask, offset_x, offset_y = extract_cutouts_vectorized_from_extension( - mock_hdu_ones, - mock_wcs, - np.array([150.0]), - np.array([2.0]), - np.array([10]), - source_ids=["flux_test"], - padding_factor=1.0, - config=config, + cutouts, success_mask, offset_x, offset_y, *_ = ( + extract_cutouts_vectorized_from_extension( + mock_hdu_ones, + mock_wcs, + np.array([150.0]), + np.array([2.0]), + np.array([10]), + source_ids=["flux_test"], + padding_factor=1.0, + config=config, + ) ) # Verify flux conversion was called @@ -210,7 +218,7 @@ def flux_conv_side_effect(config, data, header): mock_flux_conv.side_effect = flux_conv_side_effect size = 20 - cutouts, success_mask, offset_x, offset_y = ( + cutouts, success_mask, offset_x, offset_y, *_ = ( extract_cutouts_vectorized_from_extension( mock_hdu_ones, mock_wcs, @@ -239,7 +247,7 @@ def test_padding_factor_zoom_in(self, mock_hdu_gradient, mock_wcs): padding_factor = 0.5 expected_extraction_size = int(size * padding_factor) - cutouts, success_mask, offset_x, offset_y = extract_cutouts_vectorized_from_extension( + cutouts, success_mask, offset_x, offset_y, *_ = extract_cutouts_vectorized_from_extension( mock_hdu_gradient, mock_wcs, np.array([150.0]), @@ -266,7 +274,7 @@ def test_padding_factor_zoom_out(self, mock_hdu_gradient, mock_wcs): padding_factor = 2.0 expected_extraction_size = int(size * padding_factor) - cutouts, success_mask, offset_x, offset_y = extract_cutouts_vectorized_from_extension( + cutouts, success_mask, offset_x, offset_y, *_ = extract_cutouts_vectorized_from_extension( mock_hdu_gradient, mock_wcs, np.array([150.0]), @@ -294,15 +302,17 @@ def test_gradient_preservation(self, mock_hdu_gradient, mock_wcs): mock_world_to_pixel.return_value = (30.0, 40.0) size = 5 # Small size to manually verify - cutouts, success_mask, offset_x, offset_y = extract_cutouts_vectorized_from_extension( - mock_hdu_gradient, - mock_wcs, - np.array([150.0]), - np.array([2.0]), - np.array([size]), - source_ids=["gradient_test"], - padding_factor=1.0, - config=None, + cutouts, success_mask, offset_x, offset_y, *_ = ( + extract_cutouts_vectorized_from_extension( + mock_hdu_gradient, + mock_wcs, + np.array([150.0]), + np.array([2.0]), + np.array([size]), + source_ids=["gradient_test"], + padding_factor=1.0, + config=None, + ) ) assert success_mask[0], "Extraction failed" @@ -325,7 +335,7 @@ def test_batch_extraction_consistency(self, mock_hdu_ones, mock_wcs): size_array = np.array([5, 11, 20]) # Batch extraction - batch_cutouts, batch_success, batch_offset_x, batch_offset_y = ( + batch_cutouts, batch_success, batch_offset_x, batch_offset_y, *_ = ( extract_cutouts_vectorized_from_extension( mock_hdu_ones, mock_wcs, @@ -340,7 +350,7 @@ def test_batch_extraction_consistency(self, mock_hdu_ones, mock_wcs): # Individual extractions for i, size in enumerate(size_array): - single_cutouts, single_success, single_offset_x, single_offset_y = ( + single_cutouts, single_success, single_offset_x, single_offset_y, *_ = ( extract_cutouts_vectorized_from_extension( mock_hdu_ones, mock_wcs, @@ -372,7 +382,7 @@ def test_flux_conversion_bug_regression(self, mock_hdu_ones, mock_wcs): test_sizes = [5, 11, 15, 21] for size in test_sizes: - cutouts, success_mask, offset_x, offset_y = ( + cutouts, success_mask, offset_x, offset_y, *_ = ( extract_cutouts_vectorized_from_extension( mock_hdu_ones, mock_wcs, @@ -434,7 +444,7 @@ def _extract_and_verify_offset(self, hdu, wcs, target_x, target_y, cutout_size): """Helper to extract cutout and verify offset matches expected value.""" target_ra, target_dec = wcs.pixel_to_world_values(target_x, target_y) - cutouts, success_mask, offset_x, offset_y = extract_cutouts_vectorized_from_extension( + cutouts, success_mask, offset_x, offset_y, *_ = extract_cutouts_vectorized_from_extension( hdu, wcs, np.array([target_ra]), @@ -549,7 +559,7 @@ def test_batch_mixed_even_odd_sizes(self, mock_hdu, mock_wcs_precise): ra_array.append(ra) dec_array.append(dec) - cutouts, success_mask, offset_x, offset_y = extract_cutouts_vectorized_from_extension( + cutouts, success_mask, offset_x, offset_y, *_ = extract_cutouts_vectorized_from_extension( mock_hdu, mock_wcs_precise, np.array(ra_array), @@ -600,7 +610,7 @@ def test_extreme_zoom_out_requires_zero_padding(self): padding_factor = 10.0 expected_size = int(target_size * padding_factor) # 640 - cutouts, success_mask, _, _ = extract_cutouts_vectorized_from_extension( + cutouts, success_mask, _, _, *_ = extract_cutouts_vectorized_from_extension( hdu=hdu, wcs_obj=wcs_obj, ra_array=np.array([180.0]), @@ -627,7 +637,7 @@ def test_large_zoom_out_edge_padding(self): padding_factor = 5.0 expected_size = int(target_size * padding_factor) # 640 - cutouts, success_mask, _, _ = extract_cutouts_vectorized_from_extension( + cutouts, success_mask, _, _, *_ = extract_cutouts_vectorized_from_extension( hdu=hdu, wcs_obj=wcs_obj, ra_array=np.array([180.0]), @@ -651,7 +661,7 @@ def test_source_outside_image_bounds(self): """Test source completely outside image fails gracefully.""" hdu, wcs_obj = self._make_mock_hdu_wcs(512, 1000.0, 1000.0) - cutouts, success_mask, _, _ = extract_cutouts_vectorized_from_extension( + cutouts, success_mask, _, _, *_ = extract_cutouts_vectorized_from_extension( hdu=hdu, wcs_obj=wcs_obj, ra_array=np.array([180.0]), @@ -670,7 +680,7 @@ def test_very_small_cutout_sizes(self, target_size): hdu, wcs_obj = self._make_mock_hdu_wcs(512, 256.0, 256.0) hdu.data = np.random.default_rng(42).random((512, 512)).astype(np.float32) - cutouts, success_mask, _, _ = extract_cutouts_vectorized_from_extension( + cutouts, success_mask, _, _, *_ = extract_cutouts_vectorized_from_extension( hdu=hdu, wcs_obj=wcs_obj, ra_array=np.array([180.0]), @@ -697,7 +707,7 @@ def test_fractional_pixel_coordinates(self): hdu.data = data target_size = 100 - cutouts, success_mask, _, _ = extract_cutouts_vectorized_from_extension( + cutouts, success_mask, _, _, *_ = extract_cutouts_vectorized_from_extension( hdu=hdu, wcs_obj=wcs_obj, ra_array=np.array([180.0]), @@ -723,7 +733,7 @@ def test_edge_zoom_out_has_zero_padding(self): padding_factor = 2.0 expected_size = int(target_size * padding_factor) # 128 - cutouts, success_mask, _, _ = extract_cutouts_vectorized_from_extension( + cutouts, success_mask, _, _, *_ = extract_cutouts_vectorized_from_extension( hdu=hdu, wcs_obj=wcs_obj, ra_array=np.array([180.0]), diff --git a/tests/cutana/unit/test_cutout_writer_fits.py b/tests/cutana/unit/test_cutout_writer_fits.py index 6bee11a..0d4a468 100644 --- a/tests/cutana/unit/test_cutout_writer_fits.py +++ b/tests/cutana/unit/test_cutout_writer_fits.py @@ -21,11 +21,13 @@ import numpy as np import pytest +from astropy import units as u from astropy.io import fits from astropy.wcs import WCS from dotmap import DotMap from cutana import cutout_writer_fits +from cutana.cutout_extraction import extract_cutouts_vectorized_from_extension from cutana.cutout_writer_fits import ( create_wcs_header, ensure_output_directory, @@ -35,6 +37,19 @@ ) +@pytest.fixture(autouse=True) +def _clear_wcs_header_cache(): + """Isolate the module-level WCS header cache between tests. + + ``cutout_writer_fits._wcs_header_cache`` is keyed on ``id(wcs)``; across tests a + freed WCS object's id can be reused, returning a stale cached header. Clearing it + per test keeps WCS assertions deterministic regardless of execution order. + """ + cutout_writer_fits._wcs_header_cache.clear() + yield + cutout_writer_fits._wcs_header_cache.clear() + + class TestCutoutWriterFitsFunctions: """Test suite for FITS writer functions.""" @@ -84,6 +99,10 @@ def mock_cutout_data(self): "processing_timestamp": 1642678800.0, "original_tile": "euclid_tile_001.fits", "tile": "euclid_tile_001.fits", + # Extraction origin/size threaded from cutout_extraction (unresized 256 px). + "extraction_origin_x": 0, + "extraction_origin_y": 0, + "extraction_size": 256, }, } @@ -118,16 +137,26 @@ def test_generate_fits_filename(self): assert filename.endswith(".fits") def test_create_wcs_header(self, mock_cutout_data): - """Test WCS header creation.""" - wcs = mock_cutout_data["wcs_info"]["VIS"] + """CRVAL is inherited from the parent tile, not re-tangented at the source.""" + wcs = mock_cutout_data["wcs_info"]["VIS"] # parent CRVAL = [150.0, 2.0] cutout_shape = (256, 256) - header = create_wcs_header(cutout_shape, original_wcs=wcs, ra_center=150.0, dec_center=2.0) + # Deliberately offset the source from the parent CRVAL so a re-tangenting + # regression (CRVAL <- source) would be caught. The extraction origin/size are + # threaded in (as they are from cutout_extraction in the real pipeline). + header = create_wcs_header( + cutout_shape, + original_wcs=wcs, + ra_source=150.05, + dec_source=2.05, + extraction_origin_x=0, + extraction_origin_y=0, + extraction_size=256, + ) - assert "CRVAL1" in header - assert "CRVAL2" in header assert "CRPIX1" in header assert "CRPIX2" in header + # CRVAL must stay at the parent tile reference, NOT the source position. assert header["CRVAL1"] == 150.0 assert header["CRVAL2"] == 2.0 @@ -154,6 +183,71 @@ def test_write_single_fits_cutout(self, mock_cutout_data, temp_output_dir): assert "NIR-Y" in ext_names assert "NIR-H" in ext_names + def test_write_single_fits_unit_and_consvflx_headers(self, mock_cutout_data, temp_output_dir): + """UNIT and CONSVFLX from cutout_data are written to the primary header.""" + output_path = temp_output_dir / "unit_header.fits" + + cutout_data = {**mock_cutout_data, "unit": "Jy", "conserved_flux": True} + success = write_single_fits_cutout(cutout_data, str(output_path), overwrite=True) + + assert success is True + with fits.open(output_path) as hdul: + assert hdul[0].header["UNIT"] == "Jy" + assert bool(hdul[0].header["CONSVFLX"]) is True + + @pytest.mark.parametrize( + ("flux_conserved_resizing", "apply_flux_conversion", "expected_unit", "expected_consvflx"), + [ + (True, True, "Jy", True), + (True, False, "OriginalUnit", True), + (False, True, "approx Jy", False), + (False, False, "approx OriginalUnit", False), + ], + ) + def test_write_fits_batch_unit_mapping( + self, + temp_output_dir, + flux_conserved_resizing, + apply_flux_conversion, + expected_unit, + expected_consvflx, + ): + """write_fits_batch maps flux config to the UNIT/CONSVFLX headers it writes.""" + cutouts_tensor = np.random.random((1, 32, 32, 1)).astype(np.float32) + batch_data = [ + { + "cutouts": cutouts_tensor, + "metadata": [ + { + "source_id": "UnitSource_001", + "ra": 150.0, + "dec": 2.0, + "tile": "euclid_tile_001.fits", + } + ], + } + ] + + written_files = write_fits_batch( + batch_data, + str(temp_output_dir), + config=DotMap( + { + "do_only_cutout_extraction": False, + "flux_conserved_resizing": flux_conserved_resizing, + "apply_flux_conversion": apply_flux_conversion, + } + ), + file_naming_template="{source_id}_cutout.fits", + create_subdirs=False, + overwrite=True, + ) + + assert len(written_files) == 1 + with fits.open(written_files[0]) as hdul: + assert hdul[0].header["UNIT"] == expected_unit + assert bool(hdul[0].header["CONSVFLX"]) == expected_consvflx + def test_write_single_fits_with_compression(self, mock_cutout_data, temp_output_dir): """Test writing FITS with compression.""" output_path = temp_output_dir / "compressed_cutout.fits" @@ -371,17 +465,34 @@ def test_create_wcs_header_comprehensive(self): original_wcs.wcs.crpix = [50.0, 50.0] original_wcs.wcs.cdelt = [-0.0001, 0.0001] + px, py = original_wcs.world_to_pixel_values(151.0, 3.0) + x_min = int(px - 64 // 2) + y_min = int(py - 64 // 2) header = create_wcs_header( - (64, 64), original_wcs=original_wcs, ra_center=151.0, dec_center=3.0 + (64, 64), + original_wcs=original_wcs, + ra_source=151.0, + dec_source=3.0, + extraction_origin_x=int(x_min), + extraction_origin_y=int(y_min), + extraction_size=64, ) - assert header["CRPIX1"] == 32.5 # 64/2 + 0.5 (FITS 1-based center) - assert header["CRPIX2"] == 32.5 # 64/2 + 0.5 (FITS 1-based center) - assert header["CRVAL1"] == 151.0 # Updated center - assert header["CRVAL2"] == 3.0 # Updated center + # CRVAL/CTYPE are inherited from the parent tile unchanged. The projection + # is NOT re-tangented at the source position (doing so keeps the parent CD + # matrix at the wrong tangent point and rotates the cutout frame). + assert header["CRVAL1"] == 150.0 # parent CRVAL preserved + assert header["CRVAL2"] == 2.0 # parent CRVAL preserved + # CRPIX is shifted to the extraction origin so the cutout reproduces the + # parent sky mapping exactly. Verify that round-trip agreement directly. + cut_wcs = WCS(header) + for cx, cy in [(0, 0), (63, 63), (10, 50)]: + sky_parent = original_wcs.pixel_to_world_values(x_min + cx, y_min + cy) + sky_cut = cut_wcs.pixel_to_world_values(cx, cy) + assert np.allclose(sky_parent, sky_cut, atol=1e-10) # Test without original WCS but with coordinates - header = create_wcs_header((128, 128), ra_center=150.5, dec_center=2.5, pixel_scale=0.6) + header = create_wcs_header((128, 128), ra_source=150.5, dec_source=2.5, pixel_scale=0.6) assert header["WCSAXES"] == 2 assert header["CTYPE1"] == "RA---TAN" @@ -475,3 +586,200 @@ def test_error_handling_comprehensive(self, temp_output_dir): invalid_data, str(temp_output_dir / "invalid_cutout.fits"), overwrite=True ) assert success is False + + +class TestCutoutWcsFidelity: + """Regression tests for cutout WCS fidelity against the parent tile. + + Guards against the re-tangenting bug where the cutout WCS set CRVAL to the + source RA/Dec and CRPIX to the geometric centre while keeping the parent tile's + CD matrix. That rotates the cutout frame by the meridian convergence between the + tile centre and the source, giving a WCS error that GROWS with distance from the + cutout centre (order ~1" at a few-arcmin FOV for sources far from the tile centre). + + The correct construction inherits the parent CRVAL/CD and only shifts CRPIX to the + extraction origin, so every cutout pixel maps to the same sky position as the + parent tile (to numerical precision). + """ + + # Declination near the pole: meridian convergence (~ tan(Dec)) is large there, + # so the re-tangenting error the fix removes is at its most pronounced. + _TILE_DEC = -85.0 + + @staticmethod + def _euclid_like_tile_wcs(): + """A MER-VIS-like tile WCS: 19200^2, TAN, 0.1"/px, CD at centre, near the pole. + + Placed at Dec = -85 deg so meridian convergence — and hence the frame + rotation a re-tangented WCS would introduce — is strong. + """ + wcs = WCS(naxis=2) + wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"] + wcs.wcs.crval = [57.9990741, TestCutoutWcsFidelity._TILE_DEC] + wcs.wcs.crpix = [9600.0, 9600.0] # tile centre, FITS 1-based + wcs.wcs.cd = [[-2.777777777778e-05, 0.0], [0.0, 2.777777777778e-05]] + wcs.wcs.cunit = ["deg", "deg"] + wcs.pixel_shape = (19200, 19200) + return wcs + + @staticmethod + def _max_sky_error_arcsec(header, tile_wcs, x_min, y_min, resize, final_size): + """Max separation (arcsec) between cutout-WCS and parent-tile sky positions. + + Ground truth uses the cv2.resize half-pixel-centre convention: + parent_pixel = origin + (cutout_pixel + 0.5) / resize - 0.5. + """ + cut_wcs = WCS(header) + max_err = 0.0 + for frac in np.linspace(0.0, 1.0, 11): + c = frac * (final_size - 1) + parent_x = x_min + (c + 0.5) / resize - 0.5 + parent_y = y_min + (c + 0.5) / resize - 0.5 + truth = tile_wcs.pixel_to_world(parent_x, parent_y) + got = cut_wcs.pixel_to_world(c, c) + max_err = max(max_err, truth.separation(got).to(u.arcsec).value) + return max_err + + def test_far_from_centre_no_growing_error(self): + """A far-off-centre source cutout must agree with the parent tile to < 1 mas.""" + tile = self._euclid_like_tile_wcs() + # Far off-centre (~0.23 deg) but the 1800 px window stays fully on-tile + # (centre in [900, 18300]) so this isolates the rotation term, not clipping. + tx, ty = 17800.3, 17800.7 + target = tile.pixel_to_world(tx, ty) + ra_c, dec_c = target.ra.deg, target.dec.deg + + requested = 1800 # 3 arcmin at 0.1"/px + # On-tile, unresized window: origin is simply the integer window start. + px, py = tile.world_to_pixel_values(ra_c, dec_c) + x_min = int(px - requested // 2) + y_min = int(py - requested // 2) + header = create_wcs_header( + (requested, requested), + original_wcs=tile, + ra_source=ra_c, + dec_source=dec_c, + extraction_origin_x=int(x_min), + extraction_origin_y=int(y_min), + extraction_size=requested, + ) + # Parent CRVAL/CD inherited, not re-tangented at the source. + assert header["CRVAL1"] == 57.9990741 + assert header["CRVAL2"] == self._TILE_DEC + + err = self._max_sky_error_arcsec(header, tile, x_min, y_min, 1.0, requested) + assert err < 1e-3, f'cutout WCS error {err:.4f}" exceeds 1 mas tolerance' + + @pytest.mark.parametrize( + "requested,padding,final", + [ + (400, 1.0, 400), # no resize, no padding + (400, 1.5, 224), # padding > 1 + downsize + (400, 0.8, 224), # padding < 1 + downsize + (300, 2.0, 600), # padding + upsize + (128, 1.0, 128), # small, exact + ], + ) + def test_padding_and_resize_combinations(self, requested, padding, final): + """WCS stays exact across padding factors and resize ratios.""" + tile = self._euclid_like_tile_wcs() + tx, ty = 18500.3, 18500.7 + target = tile.pixel_to_world(tx, ty) + ra_c, dec_c = target.ra.deg, target.dec.deg + + # On-tile window: origin is the integer window start (no clip/pad here). + px, py = tile.world_to_pixel_values(ra_c, dec_c) + ext_size = int(requested * padding) + x_min = int(px - ext_size // 2) + y_min = int(py - ext_size // 2) + resize = final / ext_size + + header = create_wcs_header( + (final, final), + original_wcs=tile, + ra_source=ra_c, + dec_source=dec_c, + extraction_origin_x=int(x_min), + extraction_origin_y=int(y_min), + extraction_size=ext_size, + ) + + err = self._max_sky_error_arcsec(header, tile, x_min, y_min, resize, final) + assert err < 1e-3, f'padding={padding} resize={resize:.3f}: WCS error {err:.4f}" too large' + + @staticmethod + def _synthetic_tile_hdu(width, height): + """Small on-disk-style tile HDU with a high-Dec TAN WCS and known pixel_shape.""" + hdr = fits.Header() + hdr["NAXIS"] = 2 + hdr["NAXIS1"] = width + hdr["NAXIS2"] = height + hdr["CTYPE1"] = "RA---TAN" + hdr["CTYPE2"] = "DEC--TAN" + hdr["CRVAL1"] = 57.999 + hdr["CRVAL2"] = -51.5 + hdr["CRPIX1"] = width / 2.0 + hdr["CRPIX2"] = height / 2.0 + hdr["CD1_1"] = -2.7778e-05 + hdr["CD1_2"] = 0.0 + hdr["CD2_1"] = 0.0 + hdr["CD2_2"] = 2.7778e-05 + hdr["CUNIT1"] = "deg" + hdr["CUNIT2"] = "deg" + data = np.zeros((height, width), dtype=np.float32) + return fits.PrimaryHDU(data=data, header=hdr), WCS(hdr) + + @pytest.mark.parametrize( + "src_x,src_y,requested,padding", + [ + (184.0, 178.0, 60, 1.0), # window overruns the tile edge -> clipped + padded + (183.0, 177.0, 50, 1.5), # heavier clip with padding + (185.0, 178.0, 20, 1.0), # fully on-tile (no clip) + (182.0, 176.0, 50, 2.0), # clip + padding upsize + ], + ) + def test_real_extraction_marker_roundtrip(self, src_x, src_y, requested, padding): + """Drive the REAL extraction path, incl. edge clipping, and check WCS fidelity. + + A marker pixel at a known parent location is extracted through the actual + ``extract_cutouts_vectorized_from_extension`` (which clips and centre-pads at + tile edges). The cutout WCS built by ``create_wcs_header`` must map the marker's + cutout pixel back to its true sky position. This is not self-confirming: the + oracle is the parent tile WCS and the real extracted data, not the CRPIX formula. + """ + cutout_writer_fits._wcs_header_cache.clear() + width = height = 200 + hdu, tile = self._synthetic_tile_hdu(width, height) + marker_px, marker_py = 185, 178 # 0-based parent pixel + hdu.data[marker_py, marker_px] = 1000.0 + marker_sky = tile.pixel_to_world(marker_px, marker_py) + + target = tile.pixel_to_world(src_x, src_y) + ra = np.array([target.ra.deg]) + dec = np.array([target.dec.deg]) + # Capture the extraction origin the way the real pipeline threads it to the writer. + cutouts, success, _, _, origin_x, origin_y = extract_cutouts_vectorized_from_extension( + hdu, tile, ra, dec, np.array([requested], dtype=int), ["s0"], padding_factor=padding + ) + assert cutouts[0] is not None and bool(success[0]) + cut = cutouts[0] + + # The marker must have been captured for this to test WCS placement. + assert cut.max() > 0, "marker not inside extracted window; adjust test params" + + header = create_wcs_header( + cut.shape, + original_wcs=tile, + ra_source=float(ra[0]), + dec_source=float(dec[0]), + extraction_origin_x=int(origin_x[0]), + extraction_origin_y=int(origin_y[0]), + extraction_size=int(requested * padding), + ) + cut_wcs = WCS(header) + cj, ci = np.unravel_index(int(np.argmax(cut)), cut.shape) # (row=y, col=x) + got = cut_wcs.pixel_to_world(int(ci), int(cj)) + err = marker_sky.separation(got).to(u.arcsec).value + assert err < 1e-3, ( + f'real-extraction marker WCS error {err:.4f}" (clip/pad handling regressed)' + ) From 5d435a061d982d4b69e0c063ecff52b046fe42c3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pablo=20G=C3=B3mez?= Date: Wed, 8 Jul 2026 11:23:19 +0200 Subject: [PATCH 2/2] fix(pyproject): point project URLs at the public repo Homepage/Repository/Issues referenced the internal ESA-Datalabs mirror; point them at the public ESA/Cutana repository. --- pyproject.toml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 408d60d..5d72ee4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -55,9 +55,9 @@ dev = [ all = ["cutana[ui,dev]"] [project.urls] -Homepage = "https://github.com/ESA-Datalabs/cutana" -Repository = "https://github.com/ESA-Datalabs/cutana" -Issues = "https://github.com/ESA-Datalabs/cutana/issues" +Homepage = "https://github.com/ESA/Cutana" +Repository = "https://github.com/ESA/Cutana" +Issues = "https://github.com/ESA/Cutana/issues" [tool.setuptools.packages.find]