diff --git a/RMS/Astrometry/AutoPlatepar.py b/RMS/Astrometry/AutoPlatepar.py index 4a5391a43..5a4e16620 100644 --- a/RMS/Astrometry/AutoPlatepar.py +++ b/RMS/Astrometry/AutoPlatepar.py @@ -362,7 +362,7 @@ def autoFitPlatepar(dir_path, config, catalog_stars, platepar_template=None, Keyword arguments (filtering parameters): photometric_sigma: [float] Sigma threshold for photometric outlier removal (default: 2.5) - fwhm_mult: [float] Multiplier of FWHM for blend detection radius (default: 2.0) + fwhm_mult: [float] Multiplier of FWHM for blend detection radius (default: 2.8) wide_fov_search: [bool] If True, use a wide FOV search range (2° to 200°) instead of the config-based range. Used as fallback when the tight search fails. @@ -624,7 +624,7 @@ def autoFitPlatepar(dir_path, config, catalog_stars, platepar_template=None, sky_obj = CatalogStar(cat_ra, cat_dec, cat_mag) # Find closest detected star to get FWHM, SNR, saturation - fwhm, snr, saturated = 2.5, 1.0, False + fwhm, snr, saturated = 1.8, 1.0, False if len(x_data) > 0: distances = np.sqrt((x_data - img_x)**2 + (y_data - img_y)**2) closest_idx = np.argmin(distances) diff --git a/RMS/Astrometry/StarFilters.py b/RMS/Astrometry/StarFilters.py index 0f4daedfe..4abfbfadb 100644 --- a/RMS/Astrometry/StarFilters.py +++ b/RMS/Astrometry/StarFilters.py @@ -19,7 +19,10 @@ # Default filtering parameters DEFAULT_PHOTOMETRIC_SIGMA = 2.5 -DEFAULT_BLEND_FWHM_MULT = 2.0 # Multiplier of FWHM for blending detection radius +# Multiplier of FWHM for blending detection radius. The FWHM convention fix in ExtractStars +# (RMS mean of axis sigmas instead of quadrature sum) shrank reported FWHMs by sqrt(2), so this +# was rescaled 2.0 -> 2.8 to keep the same effective blend radius in pixels. +DEFAULT_BLEND_FWHM_MULT = 2.8 DEFAULT_BLEND_MAG_MARGIN = 0.3 # Margin above limiting magnitude for blend check @@ -124,7 +127,7 @@ def filterBlendedStars(paired_stars, catalog_stars, platepar, jd, lim_mag, Keyword arguments: fwhm_mult: [float] Multiplier of the star's FWHM for blend detection radius. - Default is 2.0. + Default is 2.8. mag_margin: [float] Margin above lim_mag - only consider catalog stars brighter than (lim_mag + mag_margin). Default is 0.3. verbose: [bool] Print filtering info. Default is False. diff --git a/RMS/ExtractStars.py b/RMS/ExtractStars.py index f3cb0dc92..a09943583 100644 --- a/RMS/ExtractStars.py +++ b/RMS/ExtractStars.py @@ -169,10 +169,13 @@ def extractStars(img, img_median=None, mask=None, gamma=1.0, max_star_candidates # plotStars(ff, x_arr, y_arr) - # Compute FWHM from one dimensional sigma + # Compute the FWHM from the RMS mean of the two axis sigmas, so it collapses to the standard + # 2.355*sigma for a circular star. The previous quadrature sum sqrt(sigma_x^2 + sigma_y^2) + # overestimated the FWHM by a factor of sqrt(2), and disagreed with the moment-based FWHM + # measured by SkyFit2 on manually picked stars. sigma_x_fitted = np.array(sigma_x_fitted) sigma_y_fitted = np.array(sigma_y_fitted) - sigma_fitted = np.sqrt(sigma_x_fitted**2 + sigma_y_fitted**2) + sigma_fitted = np.sqrt((sigma_x_fitted**2 + sigma_y_fitted**2)/2) fwhm = 2.355*sigma_fitted return x_arr, y_arr, amplitude, intensity, fwhm, background, snr, saturated_count diff --git a/RMS/Formats/CALSTARS.py b/RMS/Formats/CALSTARS.py index b55b37637..2d6fd6017 100644 --- a/RMS/Formats/CALSTARS.py +++ b/RMS/Formats/CALSTARS.py @@ -15,9 +15,19 @@ # along with this program. If not, see . +import math import os +# CALSTARS format version, written into the file header as "Version = N". +# - Files without a Version line (implicit version 1) carry FWHM values inflated by sqrt(2): +# they were computed as 2.355*sqrt(sigma_x^2 + sigma_y^2), the quadrature sum of the two +# axis sigmas instead of their RMS mean. +# - Version 2: the FWHM column is the standard Gaussian FWHM (2.355*sigma for a circular +# star). readCALSTARS() normalizes version 1 files to this convention on read. +CALSTARS_VERSION = 2 + + def writeCALSTARS(star_list, ff_directory, file_name, cam_code, nrows, ncols, chunk_frames=256): """ Writes the star list into the CAMS CALSTARS format. @@ -53,6 +63,7 @@ def writeCALSTARS(star_list, ff_directory, file_name, cam_code, nrows, ncols, ch star_file.write("Ncols = " + str(ncols) + "\n") star_file.write("Nframes = " + str(chunk_frames) + "\n") star_file.write("Nstars = -1" + "\n") + star_file.write("Version = {:d}".format(CALSTARS_VERSION) + "\n") # Write all stars in the CALSTARS file for star in star_list: @@ -105,12 +116,18 @@ def readCALSTARS(file_path, file_name, chunk_frames=256): Will be overwritten by a number in the CALSTARS file if present. Return: - star_list, chunk_frames: + star_list, chunk_frames: - star_list [list] a list of star data, entries: ff_name, star_data star_data entries: x, y, bg_level, level, fwhm - chunk_frames [int] Number of frames in the FF file or frame chunk. + + Note: + The FWHM values are always returned in the standard convention (2.355*sigma for a + circular star), regardless of the file version. Version 1 files (no "Version" header + line) store FWHM inflated by sqrt(2) and are normalized on read; the -1.0 sentinel + (no FWHM available) is left untouched. """ @@ -126,6 +143,10 @@ def readCALSTARS(file_path, file_name, chunk_frames=256): calibrationstars_list = [] + # Files without a Version header line are implicitly version 1 (FWHM inflated by + # sqrt(2), normalized below) + version = 1 + ff_name = '' star_data = [] skip_lines = 0 @@ -141,6 +162,11 @@ def readCALSTARS(file_path, file_name, chunk_frames=256): chunk_frames = int(line.split('=')[-1]) continue + # Read the format version if given (Version = ...) + if "Version" in line: + version = int(line.split('=')[-1]) + continue + # Check for end of star entry if (("===" in line) or ("###" in line)) and len(ff_name): @@ -181,6 +207,12 @@ def readCALSTARS(file_path, file_name, chunk_frames=256): # Read FWHM if given if len(line) >= 5: fwhm = float(line[4]) + + # Normalize version 1 FWHM values (quadrature sum of axis sigmas) to the + # standard convention, leaving the -1.0 "not available" sentinel untouched + if (version < 2) and (fwhm > 0): + fwhm = fwhm/math.sqrt(2) + else: fwhm = -1.0 diff --git a/Utils/Flux.py b/Utils/Flux.py index ece2047b9..b4eb0fa9a 100644 --- a/Utils/Flux.py +++ b/Utils/Flux.py @@ -107,8 +107,10 @@ def __init__(self): # Minimum number of meteors in the time bin self.meteors_min = 3 - # Default star FWHM, it it's not available (pz) - self.default_fwhm = 3 + # Default star FWHM if it's not available in CALSTARS, in the standard convention + # (was 3 under the old sqrt(2)-inflated FWHM convention; readCALSTARS now normalizes + # all files to the standard convention) + self.default_fwhm = 2.1 # Filter out nights which have too many detections - it is assumed that the false positives are # present if there are too many sporadic meteors @@ -1895,9 +1897,19 @@ def getSensorCharacterization(dir_path, config, flux_config, meteor_data, defaul data = " ".join(f.readlines()) sensor_data = json.loads(data) - # Remove the info entry + # Cached files without a version entry were computed under the old + # sqrt(2)-inflated FWHM convention - normalize them to the standard convention + # (readCALSTARS does the same for CALSTARS files) + if '-2' not in sensor_data: + for key in sensor_data: + if key != '-1' and sensor_data[key][0] > 0: + sensor_data[key][0] /= np.sqrt(2) + + # Remove the info entries if '-1' in sensor_data: del sensor_data['-1'] + if '-2' in sensor_data: + del sensor_data['-2'] # If file FWHM is -1 and the default FWHM is not, override it for key in sensor_data: @@ -1918,9 +1930,11 @@ def getSensorCharacterization(dir_path, config, flux_config, meteor_data, defaul # Save to file for posterior use with open(sensor_characterization_path, 'w') as f: - # Add an explanation what each entry means + # Add an explanation what each entry means, and a format version marker + # (version 2 = FWHM in the standard convention; absent = sqrt(2)-inflated) sensor_data_save = dict(sensor_data) sensor_data_save['-1'] = {"FF file name": ['median star FWHM']} + sensor_data_save['-2'] = {"version": [2]} # Convert collection areas to JSON out_str = json.dumps(sensor_data_save, indent=4, sort_keys=True) diff --git a/Utils/SkyFit2.py b/Utils/SkyFit2.py index a2b8a2aeb..d32adc4eb 100644 --- a/Utils/SkyFit2.py +++ b/Utils/SkyFit2.py @@ -179,7 +179,8 @@ def computeSolarSystemMagnitude(body_name, body, sun, time): limitingMagnitude, screenNudgeToAzAltDelta, fovCentreZenithDirection from RMS.Astrometry.AtmosphericExtinction import atmosphericExtinctionCorrection from RMS.Astrometry.StarClasses import CatalogStar, GeoPoint, PlanetPoint, PairedStars -from RMS.Astrometry.StarFilters import filterPhotometricOutliers, filterBlendedStars +from RMS.Astrometry.StarFilters import filterPhotometricOutliers, filterBlendedStars, \ + DEFAULT_BLEND_FWHM_MULT from RMS.Astrometry.Conversions import date2JD, JD2HourAngle, trueRaDec2ApparentAltAz, \ apparentAltAz2TrueRADec, J2000_JD, jd2Date, datetime2JD, JD2LST, geo2Cartesian, vector2RaDec, raDec2Vector from RMS.Astrometry.AstrometryNet import astrometryNetSolve @@ -5620,11 +5621,13 @@ def _countTrueFalsePositives(self, ff_name, intensity_threshold, segment_radius, # Count matches using KD-tree for fast nearest-neighbor lookup # Dynamic match radius based on FWHM - wide stars have less precise centroids + # (multiplier rescaled 1.5 -> 2.1 to preserve the effective radius after the + # sqrt(2) FWHM convention fix in ExtractStars) cat_tree = cKDTree(np.column_stack([catalog_x, catalog_y])) det_coords = np.column_stack([x_arr, y_arr]) nn_dist, _ = cat_tree.query(det_coords, k=1) - effective_radii = np.maximum(match_radius, 1.5 * fwhm_arr) + effective_radii = np.maximum(match_radius, 2.1 * fwhm_arr) matched = nn_dist <= effective_radii n_true_pos = int(np.sum(matched)) @@ -5749,8 +5752,10 @@ def _findSegmentRadiusFromFWHM(self, ff_name, visible_cat_x, visible_cat_y): segment_radius > 2 * sigma / sqrt(max_feature_ratio) where sigma = FWHM / 2.355. - Uses a multiplier of 1.5x the 90th percentile FWHM to provide margin for - the Gaussian fit while keeping the segment compact. + Uses a 2.1x margin on the minimum segment derived from the 90th percentile FWHM to + provide room for the Gaussian fit while keeping the segment compact (rescaled from + 1.5x to preserve the selected segment_radius after the sqrt(2) FWHM convention fix + in ExtractStars). Arguments: ff_name: [str] Name of the FF file. @@ -5813,7 +5818,9 @@ def _findSegmentRadiusFromFWHM(self, ff_name, visible_cat_x, visible_cat_y): matched_fwhms = [] n_true_pos = 0 for det_x, det_y, det_fwhm in zip(x_arr, y_arr, fwhm_arr): - effective_radius = max(match_radius, 1.5 * det_fwhm) + # Multiplier rescaled 1.5 -> 2.1 to preserve the effective radius after the + # sqrt(2) FWHM convention fix in ExtractStars + effective_radius = max(match_radius, 2.1 * det_fwhm) distances = np.sqrt((visible_cat_x - det_x)**2 + (visible_cat_y - det_y)**2) if np.min(distances) <= effective_radius: n_true_pos += 1 @@ -5846,11 +5853,12 @@ def _findSegmentRadiusFromFWHM(self, ff_name, visible_cat_x, visible_cat_y): # Compute segment_radius from FWHM # Need: segment_radius > 2 * sigma / sqrt(max_feature_ratio) # sigma = FWHM / 2.355 - # With 1.5x margin for robust fitting: + # With 2.1x margin for robust fitting (rescaled from 1.5x to preserve the selected + # segment_radius after the sqrt(2) FWHM convention fix in ExtractStars): max_feature_ratio = getattr(self.config, 'max_feature_ratio', 0.8) sigma_90 = fwhm_90 / 2.355 min_segment = 2 * sigma_90 / np.sqrt(max_feature_ratio) - best_segment = int(np.ceil(min_segment * 1.5)) + best_segment = int(np.ceil(min_segment * 2.1)) # Clamp to valid range [4, 20] best_segment = max(4, min(20, best_segment)) @@ -5858,9 +5866,9 @@ def _findSegmentRadiusFromFWHM(self, ff_name, visible_cat_x, visible_cat_y): print(f" FWHM stats ({len(all_fwhms)} stars): " f"median={fwhm_median:.1f}, 90th={fwhm_90:.1f}, max={fwhm_max:.1f} px") print(f" min segment for 90th FWHM: {min_segment:.1f} px " - f"(with 1.5x margin: {min_segment * 1.5:.1f})") + f"(with 2.1x margin: {min_segment * 2.1:.1f})") print(f"\n Selected segment_radius: {best_segment} " - f"(from 1.5x FWHM-derived minimum)") + f"(from 2.1x FWHM-derived minimum)") return best_segment, n_true_pos, fp_ratio @@ -8827,7 +8835,7 @@ def filterPositionalOutliers(self, sigma_threshold=3.0, abs_floor_px=3.0): return removed_count - def filterBlendedStars(self, fwhm_mult=2.0, mag_margin=0.3): + def filterBlendedStars(self, fwhm_mult=DEFAULT_BLEND_FWHM_MULT, mag_margin=0.3): """ Filter paired_stars by removing likely blended stars. @@ -12869,7 +12877,7 @@ def iteration_callback(iteration, pp_iter, outlier_mask, rmsd_arcmin): sky_obj = CatalogStar(cat_star[0], cat_star[1], cat_star[2]) # Look up FWHM, SNR, and saturation from detected stars - fwhm, snr, saturated = 2.5, 1.0, False + fwhm, snr, saturated = 1.8, 1.0, False if len(det_x) > 0: distances = np.sqrt((det_x - img_x)**2 + (det_y - img_y)**2) closest_idx = np.argmin(distances) @@ -12908,7 +12916,7 @@ def iteration_callback(iteration, pp_iter, outlier_mask, rmsd_arcmin): # Filter blended stars before final fit if len(self.paired_stars) >= 15: - removed = self.filterBlendedStars(fwhm_mult=2.0, mag_margin=0.3) + removed = self.filterBlendedStars(fwhm_mult=DEFAULT_BLEND_FWHM_MULT, mag_margin=0.3) if removed > 0: print("Pairs after blend filtering: {}".format(len(self.paired_stars))) @@ -13873,7 +13881,7 @@ def iteration_callback(iteration, pp_iter, outlier_mask, rmsd_arcmin): # Filter blended stars before final fit if len(self.paired_stars) >= 15: - removed = self.filterBlendedStars(fwhm_mult=2.0, mag_margin=0.3) + removed = self.filterBlendedStars(fwhm_mult=DEFAULT_BLEND_FWHM_MULT, mag_margin=0.3) if removed > 0: print("Pairs after blend filtering: {}".format(len(self.paired_stars)))