Improve fireball detector recall with stdpixel decontamination

.config

@@ -516,6 +516,15 @@ point_ratio_threshold: 0.7
 ; Maximum number of lines which are allowed to be found on the image
 max_lines: 5
 
+; Maximum angular velocity for fireball candidates in deg/s. Converted to a per-sensor
+; px/frame cap at runtime using fps and the deg/pixel scale derived from the nominal FOV
+; (fov_w/fov_h in the [Capture] section above), so the same value transfers between cameras
+; with different FOV/resolution. Keep those nominal FOV values correct for the cap to be
+; meaningful. Default is wider than the meteor detector's ang_vel_max (35 deg/s) because the
+; fireball detector is recall-first and the cap exists mostly to reject birds/planes, not to
+; constrain real fireballs.
+fireball_max_ang_vel: 60.0
+
 
 
 [MeteorDetection]

RMS/ConfigReader.py

@@ -560,6 +560,13 @@ def __init__(self):
         self.line_distance_const = 4   # constant that determines the influence of average point distance on the line quality
         self.point_ratio_threshold = 0.7# ratio of how many points must be close to the line before considering searching for another line
         self.max_lines = 5             # maximum number of lines
+        # Maximum angular velocity for fireball candidates in deg/s. Converted to a per-sensor
+        # px/frame cap at runtime using fps and the deg/pixel scale derived from the nominal FOV
+        # (fov_w/fov_h in the [Capture] section), so the same value transfers between cameras with
+        # different FOV/resolution. Default is wider than the meteor detector's ang_vel_max
+        # (35 deg/s) because the fireball detector is recall-first and the cap exists mostly to
+        # reject birds/planes, not to constrain real fireballs.
+        self.fireball_max_ang_vel = 60.0
 
         ##### MeteorDetection
 
@@ -1512,7 +1519,9 @@ def parseFireballDetection(config, parser):
 
     if parser.has_option(section, "max_lines"):
         config.max_lines = parser.getint(section, "max_lines")
-
+
+    if parser.has_option(section, "fireball_max_ang_vel"):
+        config.fireball_max_ang_vel = parser.getfloat(section, "fireball_max_ang_vel")
 
 
 

RMS/Routines/Grouping3D.py

@@ -138,52 +138,78 @@ def __init__(self, config):
 
 
 
-def findCoefficients(line_list):
+def findCoefficients(line_list, config=None):
     """ Extract coefficients from list of lines that can be consumed by RMS.VideoExtraction.
-    
+
     Arguments:
         line_list: [list] list of detected lines
-
+        config: [config object] configuration parameters (optional). Used to convert the
+            fireball_max_ang_vel limit (deg/s) into a per-sensor px/frame threshold via
+            fps and the average deg/pixel scale (fov_h/height + fov_w/width)/2.
+
     Return:
-        coeff: [list] coefficients for each detected line in format: [first point, slope of XZ, slope of YZ, 
+        coeff: [list] coefficients for each detected line in format: [first point, slope of XZ, slope of YZ,
             first frame, last frame]
     """
-
+
+    # Fallback px/frame cap matches the legacy hard-coded value when no config is available
+    # (e.g. unit tests, external callers). At default 720p/25fps/f=16 this is ~39 deg/s --
+    # tighter than the production default of 60 deg/s, so tests get conservative behavior.
+    max_velocity_px = 2.0
+
+    if config is not None:
+        max_ang_vel = getattr(config, 'fireball_max_ang_vel', None)
+        fps = getattr(config, 'fps', None)
+        fov_h = getattr(config, 'fov_h', None)
+        fov_w = getattr(config, 'fov_w', None)
+        height = getattr(config, 'height', None)
+        width = getattr(config, 'width', None)
+        f = getattr(config, 'f', None)
+
+        # Check each required attribute individually for clearer diagnostics
+        missing = [name for name, val in [('fireball_max_ang_vel', max_ang_vel),
+                                           ('fps', fps), ('fov_h', fov_h), ('fov_w', fov_w),
+                                           ('height', height), ('width', width), ('f', f)]
+                   if not val]
+
+        if not missing:
+            # deg/s -> original px/frame via the average deg/pixel scale and fps,
+            # then -> subsampled px/frame via the decimation factor f (thresholdAndSubsample
+            # divides x,y by f, so slopes are in subsampled coordinates).
+            scale = (fov_h/float(height) + fov_w/float(width))/2.0
+            if scale > 0:
+                max_velocity_px = (max_ang_vel/scale)/float(fps)/float(f)
+        else:
+            log.warning("findCoefficients: fireball_max_ang_vel conversion skipped "
+                        "(missing/zero: %s); using fallback %.2f px/frame",
+                        ', '.join(missing), max_velocity_px)
+
     coeff = []
-    
+
     for detected_line in line_list:
-        
+
         if detected_line[0][2] < detected_line[1][2]:
             point1 = np.array(detected_line[0], dtype=np.float64)
             point2 = np.array(detected_line[1], dtype=np.float64)
         elif detected_line[0][2] > detected_line[1][2]:
             point1 = np.array(detected_line[1], dtype=np.float64)
             point2 = np.array(detected_line[0], dtype=np.float64)
         else:
-        # skip if points are on the same frame (that shouldn't happen, though)
             log.debug("Points on the same frame!")
             continue
-
-        # difference between last point and first point that represent a line
+
         point3 = point2 - point1
-
-        # slope
-        slopeXZ = point3[1]/point3[2] # speed on X axis
-        slopeYZ = point3[0]/point3[2] # speed on Y axis
-
-        # length of velocity vector
+
+        slopeXZ = point3[1]/point3[2]
+        slopeYZ = point3[0]/point3[2]
+
         total = sqrt(slopeXZ**2 + slopeYZ**2)
-
-        #print('Fireball slope:', total)
 
-        # ignore line if too fast
-        # TODO: this limit should be read from config file and calculated for FOV
-        # 1.6 is better estimate on upper speed limit, set to 2 for safety
-        if total > 2:
+        if total > max_velocity_px:
             continue
-        
-        coeff.append([point1, slopeXZ, slopeYZ, detected_line[4], detected_line[5]]) #first point, slope of XZ, slope of YZ, first frame, last frame
-        
+
+        coeff.append([point1, slopeXZ, slopeYZ, detected_line[4], detected_line[5]])
+
     return coeff
 
 

RMS/Routines/Grouping3Dcy.pyx

@@ -487,9 +487,14 @@ def thresholdAndSubsample(np.ndarray[UINT8_TYPE_t, ndim=3] frames, \
             # Compute the threshold limit
             avg_std = int(float(compressed[2, y, x]) + k1*float(compressed[3, y, x])) + j1
 
-            # Make sure the threshold limit is not above the maximum possible value
-            if avg_std > 255:
-                avg_std = 255
+            # The threshold can exceed the uint8 range and must be clipped. With the old
+            # clip at 255 a saturated pixel (max_val == 255) still passed (255 >= 255) --
+            # the real recall loss was unsaturated trail pixels whose threshold clipped to
+            # 255, and that is fixed upstream by stdpixel decontamination in
+            # VideoExtraction.findPoints. Clipping to 254 additionally lets max_val == 254
+            # pixels through: a minor, harmless extra rescue.
+            if avg_std > 254:
+                avg_std = 254
 
             if ((max_val > min_level) and (max_val >= avg_std)):
 

RMS/VideoExtraction.py

@@ -16,6 +16,8 @@
 
 
 
+from __future__ import division
+
 import math
 import time
 from multiprocessing import Process, Event
@@ -33,6 +35,93 @@
 log = getLogger("rmslogger")
 
 
+# --- Fireball stdpixel decontamination tuning constants ---
+# A bright fireball inflates stdpixel along its trail within the 256-frame FF block.
+# The fireball threshold (avepixel + k1*stdpixel + j1) then exceeds the uint8 range and
+# gets clipped, so thresholdAndSubsample keeps only the near-saturated trail core and
+# drops the fainter trail edges -- exactly where recall is lost. These constants control
+# how that self-contamination is detected and neutralised before thresholding.
+BRIGHT_PERCENTILE = 90            # maxpixel percentile above which a pixel is a bright/trail candidate
+CONTAMINATION_STD_FACTOR = 3.0    # stdpixel above this * background median counts as contaminated
+MIN_BACKGROUND_PIXELS = 100       # need at least this many background pixels to trust their median
+
+
+def decontaminateStdpixel(compressed, mask_img=None):
+    """ Replace stdpixel at fireball-contaminated pixels with the background median.
+
+    A bright fireball inflates stdpixel along its trail, pushing the detection threshold
+    (avepixel + k1*stdpixel + j1) past the uint8 ceiling where it is clipped.
+    thresholdAndSubsample then keeps only the near-saturated trail core and drops the
+    fainter edges. Replacing the inflated stdpixel with the frame's background median lets
+    the whole trail threshold normally, which is the change that actually recovers recall.
+
+    Pure numpy, no Cython, Python 2/Pi compatible. The 4-plane compressed array is copied
+    only when contamination is actually present, so the common (no-fireball) path is
+    zero-copy.
+
+    Arguments:
+        compressed: [ndarray] (4, H, W) FTP-compressed frames
+            (maxpixel, maxframe, avepixel, stdpixel).
+
+    Keyword arguments:
+        mask_img: [ndarray or None] mask image (H, W); pixels where it is 0 are excluded
+            from the background and contamination sets. Ignored when its shape does not
+            match the compressed frames (e.g. a binned mask on a full-resolution
+            extractor), mirroring maskImage()'s silent shape guard.
+
+    Return:
+        [ndarray] compressed array with stdpixel decontaminated. The input array is
+            returned unchanged (same object) when no contamination is detected.
+    """
+
+    maxpix = compressed[0]
+    stdpix = compressed[3]
+
+    # Only trust the mask when it matches the frame resolution. Extractor.__init__ bins
+    # the mask via binImageCalibration(), but the extractor runs on capture-resolution
+    # frames, so on stations with detection_binning_factor > 1 the shapes disagree. Skip
+    # the mask then instead of crashing on the boolean-and below (mirrors maskImage()).
+    if (mask_img is not None) and (mask_img.shape == maxpix.shape):
+        mask_valid = mask_img > 0
+    else:
+        mask_valid = None
+
+    # Brightness cut for the background sample. executeAll() fills masked-out pixels of
+    # maxpixel with the image mean before this runs, so restrict the percentile to
+    # mask-valid pixels to avoid skewing it on heavily-masked cameras.
+    maxpix_valid = maxpix[mask_valid] if mask_valid is not None else maxpix.ravel()
+    if maxpix_valid.size == 0:
+        return compressed
+
+    # O(n) partition instead of an O(n log n) full sort for the percentile
+    n_bright = max(1, (maxpix_valid.size*(100 - BRIGHT_PERCENTILE))//100)
+    bright_threshold = np.partition(maxpix_valid, -n_bright)[-n_bright]
+
+    bright_mask = maxpix >= bright_threshold
+    bg_mask = ~bright_mask
+    if mask_valid is not None:
+        bg_mask &= mask_valid
+
+    if np.count_nonzero(bg_mask) > MIN_BACKGROUND_PIXELS:
+        std_ref = float(np.median(stdpix[bg_mask]))
+    else:
+        std_ref = float(np.median(stdpix))
+
+    contaminated = bright_mask & (stdpix > CONTAMINATION_STD_FACTOR*max(std_ref, 1))
+    if mask_valid is not None:
+        contaminated &= mask_valid
+
+    if not np.any(contaminated):
+        return compressed
+
+    # np.array() deep-copies all 4 planes; only plane 3 (stdpixel) is then modified.
+    compressed = np.array(compressed)
+    replacement = max(1, int(round(min(std_ref, 255))))
+    compressed[3][contaminated] = np.uint8(replacement)
+
+    return compressed
+
+
 class Extractor(Process):
     """ Detects fireballs and brighter meteors on the FF files, and extracts raw frames while they are still
         in memory. 
@@ -68,8 +157,13 @@ def findPoints(self):
             (y, x, z): [tuple] Coordinates of points that form the fireball, where Z is the frame number.
         """
 
+        # Decontaminate stdpixel before thresholding so bright fireball trails are not
+        # suppressed by their own inflated stdpixel (see decontaminateStdpixel).
+        mask_img = self.mask.img if self.mask is not None else None
+        compressed = decontaminateStdpixel(self.compressed, mask_img)
+
         # Threshold and subsample frames
-        length, x, y, z = Grouping3D.thresholdAndSubsample(self.frames, self.compressed, \
+        length, x, y, z = Grouping3D.thresholdAndSubsample(self.frames, compressed, \
             self.config.min_level, self.config.min_pixels, self.config.k1, self.config.j1, self.config.f)
 
 
@@ -302,7 +396,7 @@ def executeAll(self):
         t = time.time()
 
         # Find the parameters of the line in the 3D point cloud
-        coeff = Grouping3D.findCoefficients(line_list)
+        coeff = Grouping3D.findCoefficients(line_list, self.config)
         log.debug("[" + self.filename + "] Time for finding coefficients: " + str(time.time() - t) + "s")
 
         if len(coeff) == 0: