Fix star FWHM convention: collapse to standard 2.355*sigma for circular stars

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)

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.

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

RMS/Formats/CALSTARS.py

@@ -15,9 +15,19 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 
+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
 

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)

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,21 +5853,22 @@ 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))
 
             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)))