Add the side band pedestal and the ring mean as background-corrected sums

README.md

@@ -120,6 +120,37 @@ defines the LED search areas, the first image of the experiment
 (`first_img_experiment_id`), which defines the experiment start time, and the time
 reference image.
 
+### Background correction (optional)
+
+By default the extinction coefficients are computed from `sum_col_val`, the plain sum of
+the pixel values in the search area of each LED. Step 3 additionally writes
+background-corrected sums, which can be selected with `reference_property` in
+`config_analysis.ini`:
+
+| `reference_property` | Background estimate |
+|---|---|
+| `sum_col_val` (default) | none |
+| `bgsub_sum_col_val` | median of the pixels > 0 in a 2 px border ring of the search area, subtracted pixel-wise |
+| `bgsub2_sum_col_val` | mean over the channel pixels (including 0) of the same border ring, subtracted as N·b |
+| `sidebg_sum_col_val` | mean over the channel pixels in side bands next to the LED strip, subtracted as N·B |
+
+N is the number of channel pixels in the search area. For small search areas the border
+ring lies in the halo of the LED itself, whereas the side bands sample the additive
+pedestal next to the LED: stray light, light from the smoke and noise clipped at 0. Under
+dense smoke the box sum consists mostly of this pedestal.
+
+The side bands run parallel to the strip, 25 to 60 px across it (`sideband_distances` in
+`config.ini`); pixels close to other search areas are excluded. The direction of the strip
+is estimated for each LED from its nearest neighbours in the same LED array.
+`sidebg_pedestal` holds the subtracted N·B. `sidebg_sum_col_val` is NaN if the direction
+cannot be estimated (an LED array with a single LED) or if the side bands contain no pixel,
+for example because the image ends closer to the strip. The analysis stops with an error
+naming these LEDs if the selected `reference_property` contains NaN.
+
+`bgsub2_sum_col_val` and `sidebg_sum_col_val` are kept at least at `sideband_floor`
+(default 5 %) of `sum_col_val`; values at this limit carry no information on the
+extinction.
+
 ## Demo
 
 LEDSA includes a demo that demonstrates its functionality using sample data.

doc/cli_documentation.rst

@@ -34,7 +34,7 @@ Data Extraction
      - Step 3: Analyse all images to extract pixel values inside the search areas and fit them to an LED model (SLOW!).
      - ``--r`` (Use the red channel, default), ``--g`` (Use the green channel), ``--b`` (Use the blue channel), ``-rgb`` (Run for each channel)
    * - ``-s3_fast``, ``--step_3_fast`, ``--analyse_photo_fast``
-     - Step 3: Analyse all images to extract pixel values and calculate the accumulated values inside the search areas (FAST!).
+     - Step 3: Analyse all images to extract pixel values and calculate the accumulated values inside the search areas (FAST!). Besides the plain sum, the background-corrected sums ``bgsub_sum_col_val``, ``bgsub2_sum_col_val`` and ``sidebg_sum_col_val`` are written (see ``reference_property`` in the analysis configuration).
      - --
    * - ``-re``, ``--restart``
      - Restart step 3 if it was previously interrupted. Only the images that have not been analysed are taken into account.

doc/data_extraction.rst

@@ -7,6 +7,14 @@ Data Extraction
    :show-inheritance:
 
 .. automodule:: ledsa.data_extraction.LEDAnalysisData
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+Side Band Pedestal
+------------------
+
+.. automodule:: ledsa.data_extraction.sideband_background
    :members:
    :undoc-members:
    :show-inheritance:

ledsa/analysis/ExtinctionCoefficients.py

@@ -128,6 +128,7 @@ def load_img_data(self) -> None:
         """
         Load processed image data from binary file
 
+        :raises ValueError: If the reference property contains NaN values.
         """
         if self.average_images:
             img_data = read_hdf_avg(self.experiment.channel, path=self.experiment.path)
@@ -138,6 +139,16 @@ def load_img_data(self) -> None:
         self.calculated_img_data = img_data_cropped[img_data_cropped['led_array_id'] == self.experiment.led_array]
         if self.calculated_img_data.empty:
             exit(f"Apparently there are no intensity values for led array {self.experiment.led_array}!")
+        # The solvers cannot handle NaN: the linear one stops with an unspecific error, the nonlinear one returns
+        # zeros without a warning. sidebg_sum_col_val is NaN for LEDs without side band pixels, for example.
+        nan_values = self.calculated_img_data[self.reference_property].isna()
+        if nan_values.any():
+            nan_led_ids = sorted(self.calculated_img_data[nan_values].index.get_level_values('led_id').unique())
+            led_ids_text = ', '.join(map(str, nan_led_ids[:10])) + (', ...' if len(nan_led_ids) > 10 else '')
+            raise ValueError(f'{int(nan_values.sum())} values of {self.reference_property} of led array '
+                             f'{self.experiment.led_array} are NaN (led ids {led_ids_text}). For sidebg_sum_col_val '
+                             f'this means that no side band pixels were found for these LEDs, see sideband_distances '
+                             f'in config.ini.')
         # Images without analysis results are missing from the binary data, so the IDs are not necessarily a
         # gapless range. They are kept to pair every computed coefficient with the time of its own image.
         self.analysed_img_ids = sorted(

ledsa/core/ConfigData.py

@@ -17,7 +17,8 @@ def __init__(self, load_config_file=True, img_directory=None, search_area_radius
                  time_diff_to_image_time=None, img_name_string=None, num_img_overflow=None,
                  first_img_experiment_id=None, last_img_experiment_id=None, ref_img_id=None, ignore_led_indices=None,
                  led_array_edge_indices=None, led_array_edge_coordinates=None, first_img_analysis_id=None,
-                 last_img_analysis_id=None, num_skip_imgs=0, num_skip_leds=0, merge_led_array_indices=None):  # TODO: merge LED arrays
+                 last_img_analysis_id=None, num_skip_imgs=0, num_skip_leds=0, merge_led_array_indices=None,
+                 sideband_distances='25 60', sideband_floor=0.05, sideband_exclusion_margin=None):  # TODO: merge LED arrays
         """
         :param load_config_file: Determines whether to load the config file on initialization. Defaults to True.
         :type load_config_file: bool
@@ -76,6 +77,15 @@ def __init__(self, load_config_file=True, img_directory=None, search_area_radius
         :type num_skip_leds: int
         :param merge_led_array_indices: Flag to merge LED arrays for analysis. Defaults to None. TODO: not a flag but list of arrays to merge?
         :type merge_led_array_indices: bool or None
+        :param sideband_distances: Inner and outer distance in pixels across the LED strip of the side bands used to
+            estimate the pedestal for sidebg_sum_col_val. Defaults to '25 60'.
+        :type sideband_distances: str
+        :param sideband_floor: Lower limit of sidebg_sum_col_val and bgsub2_sum_col_val as fraction of sum_col_val.
+            Defaults to 0.05.
+        :type sideband_floor: float
+        :param sideband_exclusion_margin: Side band pixels within search_area_radius + margin of any other search area
+            are excluded. Defaults to None, i.e. margin = search_area_radius.
+        :type sideband_exclusion_margin: int or None
         """
         cp.ConfigParser.__init__(self, allow_no_value=True)
         if load_config_file:
@@ -150,6 +160,13 @@ def __init__(self, load_config_file=True, img_directory=None, search_area_radius
             self['analyse_photo']['   num_skip_imgs'] = str(num_skip_imgs)
             self.set('analyse_photo', '   # In the analysis, only LEDs are considered with id dividable by skip_leds + 1.')
             self['analyse_photo']['   num_skip_leds'] = str(num_skip_leds)
+            self.set('analyse_photo', '   # Pedestal estimation in side bands parallel to the LED strip (sidebg_sum_col_val):')
+            self.set('analyse_photo', '   # inner and outer distance of the side bands across the strip in pixels,')
+            self['analyse_photo']['   sideband_distances'] = str(sideband_distances)
+            self.set('analyse_photo', '   # lower limit of the corrected sum as fraction of the plain sum,')
+            self['analyse_photo']['   sideband_floor'] = str(sideband_floor)
+            self.set('analyse_photo', '   # pixels within search_area_radius + margin of other search areas are excluded, None: margin = radius.')
+            self['analyse_photo']['   sideband_exclusion_margin'] = str(sideband_exclusion_margin)
 
             with open('config.ini', 'w', encoding='utf-8') as configfile:
                 self.write(configfile)

ledsa/core/file_handling.py

@@ -289,15 +289,19 @@ def create_binary_data(channel: int, path: str = '.') -> None:
 # prepended when they are loaded, so a file with n columns has n + 1 names.
 BASE_COLUMNS = ["img_id", "led_id", "led_array_id", "sum_col_val", "mean_col_val", "max_col_val"]
 BGSUB_COLUMNS = ["bgsub_sum_col_val"]
+SIDEBG_COLUMNS = ["sidebg_sum_col_val", "sidebg_pedestal", "bgsub2_sum_col_val"]
 LED_CENTER_COLUMNS = ["led_center_x", "led_center_y"]
 FIT_COLUMNS = ["x", "y", "dx", "dy", "A", "alpha", "wx", "wy", "fit_success", "fit_fun", "fit_nfev", "fit_time"]
 
 # known layouts by the number of names (including img_id)
 COLUMN_LAYOUTS = {
     len(BASE_COLUMNS): BASE_COLUMNS,  # -s3_fast, before the local background subtraction
-    len(BASE_COLUMNS + BGSUB_COLUMNS): BASE_COLUMNS + BGSUB_COLUMNS,  # -s3_fast
+    len(BASE_COLUMNS + BGSUB_COLUMNS): BASE_COLUMNS + BGSUB_COLUMNS,  # -s3_fast, before the side band pedestal
+    len(BASE_COLUMNS + BGSUB_COLUMNS + SIDEBG_COLUMNS): BASE_COLUMNS + BGSUB_COLUMNS + SIDEBG_COLUMNS,  # -s3_fast
     len(BASE_COLUMNS + BGSUB_COLUMNS + LED_CENTER_COLUMNS + FIT_COLUMNS):
-        BASE_COLUMNS + BGSUB_COLUMNS + LED_CENTER_COLUMNS + FIT_COLUMNS,  # -s3 with fit
+        BASE_COLUMNS + BGSUB_COLUMNS + LED_CENTER_COLUMNS + FIT_COLUMNS,  # -s3 with fit, before the side band pedestal
+    len(BASE_COLUMNS + BGSUB_COLUMNS + SIDEBG_COLUMNS + LED_CENTER_COLUMNS + FIT_COLUMNS):
+        BASE_COLUMNS + BGSUB_COLUMNS + SIDEBG_COLUMNS + LED_CENTER_COLUMNS + FIT_COLUMNS,  # -s3 with fit
 }
 # Fitted files written before the delimiter between the main and the fit part was added: max_col_value and
 # led_center_x are glued into one field, so the file has one column less and the names are shifted from there on.

ledsa/core/image_reading.py

@@ -38,6 +38,53 @@ def read_channel_data_from_img(filename: str, channel: int) -> np.ndarray:
     reader = _get_reader(filename, _read_channel_data_from_img_file, _read_channel_data_from_raw_file)
     return _read_or_raise(filename, reader, channel)
 
+
+def read_channel_data_and_mask_from_img(filename: str, channel: int) -> Tuple[np.ndarray, np.ndarray]:
+    """
+    Returns the same 2D channel array as read_channel_data_from_img together with a boolean mask of the
+    pixels that belong to the selected color channel.
+
+    For RAW files the mask is taken from the Bayer filter pattern (raw_colors_visible) and not from the pixel
+    values, because pixels of the selected channel can be 0 as well (noise clipped at 0 after the black level
+    subtraction). For JPG/PNG files every pixel belongs to the channel.
+
+    :param filename: The path of the image file to read.
+    :type filename: str
+    :param channel: The color channel to process.
+    :type channel: int
+    :return: The channel array and the boolean channel mask of the same shape.
+    :rtype: tuple[np.ndarray, np.ndarray]
+    :raises CorruptedImageError: If the image cannot be read or has an unsupported file extension.
+    """
+    def read_img_file(img_filename, img_channel):
+        channel_array = _read_channel_data_from_img_file(img_filename, img_channel)
+        return channel_array, np.ones(channel_array.shape, dtype=bool)
+
+    def read_raw_file(raw_filename, raw_channel):
+        img_array, filter_array = _read_img_array_from_raw_file(raw_filename, raw_channel)
+        channel_mask = get_bayer_channel_mask(filter_array, raw_channel)
+        return np.where(channel_mask, img_array, 0), channel_mask
+
+    reader = _get_reader(filename, read_img_file, read_raw_file)
+    return _read_or_raise(filename, reader, channel)
+
+
+def get_bayer_channel_mask(filter_array: np.ndarray, channel: int) -> np.ndarray:
+    """
+    Returns a boolean mask of the Bayer filter positions that belong to the given color channel.
+    The Bayer array contains two green channels (1 and 3), one red (0) and one blue (2).
+
+    :param filter_array: Bayer filter pattern of the raw image (rawpy raw_colors_visible).
+    :type filter_array: np.ndarray
+    :param channel: The color channel (0 red, 1 green, 2 blue).
+    :type channel: int
+    :return: Boolean mask, True where the pixel belongs to the channel.
+    :rtype: np.ndarray
+    """
+    if channel == 1:
+        return (filter_array == 1) | (filter_array == 3)
+    return filter_array == channel
+
 def read_img_array_from_img(filename: str, channel: int) -> np.ndarray:
     """
     Returns a 2D array of the image for all color channels.

ledsa/data_extraction/LEDAnalysisData.py

@@ -21,6 +21,16 @@ class LEDAnalysisData:
     :ivar bgsub_sum_color_value: Integrated color value of the LED over the search area after subtraction of the
         local background estimated from the border of the search area.
     :vartype bgsub_sum_color_value: float
+    :ivar sidebg_sum_color_value: Integrated color value of the LED over the search area after subtraction of the
+        pedestal estimated in side bands next to the LED strip, at least a fraction of the plain sum.
+    :vartype sidebg_sum_color_value: float
+    :ivar sidebg_pedestal: Pedestal N * B subtracted for sidebg_sum_color_value (N channel pixels of the search
+        area, B mean per channel pixel in the side bands).
+    :vartype sidebg_pedestal: float
+    :ivar bgsub2_sum_color_value: Integrated color value of the LED over the search area minus N * b, with b the mean
+        over the channel pixels (including 0) of the border ring of the search area, at least a fraction of the
+        plain sum (sideband_floor).
+    :vartype bgsub2_sum_color_value: float
     :ivar fit_results: Fit results after fitting.
     :vartype fit_results: OptimizeResult
     :ivar fit_time: Time taken for fitting.
@@ -45,6 +55,9 @@ def __init__(self, led_id, led_array, fit_leds):
         self.sum_color_value = None
         self.max_color_value = None
         self.bgsub_sum_color_value = None
+        self.sidebg_sum_color_value = None
+        self.sidebg_pedestal = None
+        self.bgsub2_sum_color_value = None
         self.fit_leds = fit_leds
         self.fit_results = None
         self.fit_time = None
@@ -73,6 +86,8 @@ def get_main_data_string(self) -> str:
         out_str = f'{self.led_id:4d},{self.led_array:2d},'
         out_str += f'{self.sum_color_value:10.4e},{self.mean_color_value:10.4e},{self.max_color_value}'
         out_str += f',{self.bgsub_sum_color_value:10.4e}'
+        out_str += f',{self.sidebg_sum_color_value:10.4e},{self.sidebg_pedestal:10.4e}'
+        out_str += f',{self.bgsub2_sum_color_value:10.4e}'
         return out_str
 
     def get_fit_data_string(self) -> str: