Feature/local background subtraction

ledsa/analysis/ConfigDataAnalysis.py

@@ -31,6 +31,8 @@ def __init__(self, load_config_file=True, camera_position=None, num_layers=20, d
         :param num_cores: Number of CPU cores for (multicore) processing. If greater than 1, multicore processing is applied. Defaults to 1.
         :type num_cores: int
         :param reference_property: Property used for reference in LEDSA. Defaults to 'sum_col_val'.
+            Use 'bgsub_sum_col_val' to base the analysis on intensities with the local background
+            subtracted, estimated from the border of each LED search area.
         :type reference_property: str
         :param average_images: Determines if intensities are computed as an average from two consecutive images. Defaults to False.
         :type average_images: bool

ledsa/core/file_handling.py

@@ -293,6 +293,8 @@ def _get_column_names(channel: int) -> List[str]:
     parameters = ledsa.core.file_handling.read_table(file_path, delim=',', silent=True)
     columns = ["img_id", "led_id", "led_array_id",
                "sum_col_val", "mean_col_val", "max_col_val"]
+    if parameters.shape[1] in (6, 20):  # files written since the local background subtraction was added
+        columns.append("bgsub_sum_col_val")
     if parameters.shape[1] > len(columns):
         columns.extend(["led_center_x", "led_center_y"])
         columns.extend(["x", "y", "dx", "dy", "A", "alpha", "wx", "wy", "fit_success", "fit_fun", "fit_nfev"])

ledsa/data_extraction/LEDAnalysisData.py

@@ -18,6 +18,9 @@ class LEDAnalysisData:
     :vartype sum_color_value: float
     :ivar max_color_value: Maximum color value observed for the LED.
     :vartype max_color_value: float
+    :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 fit_results: Fit results after fitting.
     :vartype fit_results: OptimizeResult
     :ivar fit_time: Time taken for fitting.
@@ -41,6 +44,7 @@ def __init__(self, led_id, led_array, fit_leds):
         self.mean_color_value = None
         self.sum_color_value = None
         self.max_color_value = None
+        self.bgsub_sum_color_value = None
         self.fit_leds = fit_leds
         self.fit_results = None
         self.fit_time = None
@@ -67,6 +71,7 @@ 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}'
         return out_str
 
     def get_fit_data_string(self) -> str:

ledsa/data_extraction/step_3_functions.py

@@ -158,10 +158,40 @@ def _generate_led_analysis_data(conf: ConfigData, channel: int, data: np.ndarray
     led_data.mean_color_value = np.mean(data[search_area])
     led_data.sum_color_value = np.sum(data[search_area])
     led_data.max_color_value = np.amax(data[search_area])
+    background = _estimate_local_background(data[search_area])
+    led_data.bgsub_sum_color_value = float(np.sum(np.clip(data[search_area].astype(np.float64) - background,
+                                                          0.0, None)))
 
     return led_data
 
 
+def _estimate_local_background(search_area_data: np.ndarray, border_width: int = 2) -> float:
+    """
+    Estimate the local background of a search area from the median of its border pixels.
+
+    The LED sits at the center of the search area, so the outer border ring samples the
+    scene background (stray light, smoke path radiance). Pixels with value 0 are excluded
+    because raw Bayer arrays mask all pixels of foreign color channels with 0 and the
+    black level subtraction clips at 0.
+
+    :param search_area_data: Part of the image where the LED is located.
+    :type search_area_data: np.ndarray
+    :param border_width: Width of the border ring in pixels.
+    :type border_width: int
+    :return: Estimated background value per pixel. 0 if no valid border pixels exist.
+    :rtype: float
+    """
+    if search_area_data.shape[0] <= 2 * border_width or search_area_data.shape[1] <= 2 * border_width:
+        return 0.0
+    border_mask = np.ones(search_area_data.shape, dtype=bool)
+    border_mask[border_width:-border_width, border_width:-border_width] = False
+    border_pixels = search_area_data[border_mask]
+    border_pixels = border_pixels[border_pixels > 0]
+    if border_pixels.size == 0:
+        return 0.0
+    return float(np.median(border_pixels))
+
+
 def _save_results_in_file(channel: int, img_data: LEDAnalysisData, img_filename: str, img_id: str, img_infos: np.ndarray, basename: str) -> None:
     """
     Save analysis results to a file.
@@ -302,7 +332,7 @@ def _create_header(channel: int, img_id: str, img_filename: str, img_infos: np.n
     out_str = f'# image root = {basename}, photo file name = {img_filename}, '
     out_str += f"channel = {channel}, "
     out_str += f"time[s] = {img_infos[int(img_id) - 1][3]}\n"
-    out_str += "# id,line,sum_col_value,average_col_value,max_col_value"
+    out_str += "# id,line,sum_col_value,average_col_value,max_col_value,bgsub_sum_col_value"
     if fit_leds:
         out_str += ",led_center_x, led_center_y"
         out_str += ",x,y,dx,dy,A,alpha,wx,wy,fit_success,fit_fun,fit_nfev,fit_time"

ledsa/tests/AcceptanceTests/07_test_background_subtraction.robot

@@ -0,0 +1,70 @@
+*** Settings ***
+Resource  global_keywords.resource
+
+Force Tags  analysis  background_subtraction
+
+*** Variables ***
+${EXTINCTION_DIR}    analysis${/}extinction_coefficients${/}linear
+${SUM_RESULT}        ${EXTINCTION_DIR}${/}extinction_coefficients_linear_channel_0_sum_col_val_led_array_0.csv
+${SUM_REFERENCE}     ${EXTINCTION_DIR}${/}extinction_coefficients_linear_channel_0_sum_col_val_led_array_0_before_strip.csv
+
+*** Test Cases ***
+Bgsub Reference Property Matches Ground Truth On Clean Images
+    [Documentation]    On background-free test images there is nothing to subtract, so the
+    ...                opt-in bgsub property must reproduce the ground truth just like the
+    ...                default property does.
+    Change Directory    ${WORKDIR}
+    Create And Fill Config Analysis    linear    bgsub_sum_col_val
+    Execute Ledsa    --analysis
+    Check Bgsub Results    2
+    Check Bgsub Results    3
+    Check Bgsub Results    4
+
+Default Analysis Is Not Changed By The Bgsub Column
+    [Documentation]    Results with the default reference property must be identical whether
+    ...                the step 3 files contain the new bgsub column or not (old data format).
+    Change Directory    ${WORKDIR}
+    Create And Fill Config Analysis    linear
+    Execute Ledsa    --analysis
+    Copy File    ${SUM_RESULT}    ${SUM_REFERENCE}
+    Strip Bgsub Column From Led Position Files
+    Execute Ledsa    --analysis
+    ${identical} =    Check Extinction Coefficient Files Are Identical    ${SUM_RESULT}    ${SUM_REFERENCE}
+    Should Be True    ${identical}
+
+Background Light Biases Sum But Not Bgsub
+    [Documentation]    With a constant additive background in the images the default property
+    ...                underestimates the extinction while the bgsub property stays close to
+    ...                the ground truth.
+    Create Directory    ${WORKDIR}${/}background
+    Change Directory    ${WORKDIR}${/}background
+    Create Test Data    background=15
+    Create Config
+    Execute Ledsa    -s1
+    Execute Ledsa    -s2
+    Execute Ledsa    --coordinates
+    Execute Ledsa    -s3_fast
+    Create And Fill Config Analysis    linear
+    Execute Ledsa    --analysis
+    Create And Fill Config Analysis    linear    bgsub_sum_col_val
+    Execute Ledsa    --analysis
+    Check Background Results    2
+    Check Background Results    3
+    Check Background Results    4
+
+*** Keywords ***
+Check Bgsub Results
+    [Arguments]    ${image_id}
+    ${rmse} =    Check Input Vs Computed Extinction Coefficients    ${image_id}    linear
+    ...    reference_property=bgsub_sum_col_val
+    # slightly more slack than the 0.05 of the default property: the border ring of the
+    # synthetic JPGs contains faint LED tails and compression artefacts that get subtracted
+    Should Be True    ${rmse} < 0.06
+
+Check Background Results
+    [Arguments]    ${image_id}
+    ${rmse_sum} =    Check Input Vs Computed Extinction Coefficients    ${image_id}    linear
+    ${rmse_bgsub} =    Check Input Vs Computed Extinction Coefficients    ${image_id}    linear
+    ...    reference_property=bgsub_sum_col_val
+    Should Be True    ${rmse_sum} > 0.05
+    Should Be True    ${rmse_bgsub} < 0.05

ledsa/tests/AcceptanceTests/LedsaATestLibrary.py

@@ -1,3 +1,4 @@
+import glob
 import os
 from subprocess import Popen, PIPE
 
@@ -24,7 +25,7 @@ def change_dir(self, new_dir):
         os.chdir(new_dir)
 
     @keyword
-    def create_test_data(self, num_of_leds=100, num_of_layers=20, bottom_border=0, top_border=3):
+    def create_test_data(self, num_of_leds=100, num_of_layers=20, bottom_border=0, top_border=3, background=0):
         # Create test_data directory if it doesn't exist
         if not os.path.exists('test_data'):
             os.makedirs('test_data')
@@ -65,7 +66,7 @@ def extco_quad(z):
             for z in np.linspace(bottom_border + 0.05, top_border - 0.05, num_of_leds):
                 ex.add_led(0, 4, z)
             ex.set_extinction_coefficients(extinction_coefficients)
-            create_test_image(image_id, ex)
+            create_test_image(image_id, ex, background=int(background))
 
     # ------------------------------------------------------------------
     # Stacked / multi-camera keywords
@@ -270,8 +271,10 @@ def plot_input_vs_computed_extinction_coefficients(self, solver, first=1, last=4
             plt.close()
 
     @keyword
-    def check_input_vs_computed_extinction_coefficients(self, image_id, solver, led_array=0, channel=0):
-        _, extinction_coefficients_computed = load_extinction_coefficients_computed(solver, channel, led_array)
+    def check_input_vs_computed_extinction_coefficients(self, image_id, solver, led_array=0, channel=0,
+                                                        reference_property='sum_col_val'):
+        _, extinction_coefficients_computed = load_extinction_coefficients_computed(solver, channel, led_array,
+                                                                                    reference_property)
         extinction_coefficients_input = np.loadtxt(os.path.join('test_data', f'test_extinction_coefficients_input_{image_id}.csv'), delimiter=',')
         rmse = np.sqrt(
             np.mean((extinction_coefficients_input - extinction_coefficients_computed[int(image_id) - 1, :]) ** 2))
@@ -297,10 +300,10 @@ def create_and_fill_config(self, first=1, last=4):
         conf.save()
 
     @keyword
-    def create_and_fill_config_analysis(self, solver):
+    def create_and_fill_config_analysis(self, solver, reference_property='sum_col_val'):
         conf = ConfigDataAnalysis(load_config_file=False, camera_position=None, num_layers=20, domain_bounds=None,
                                   led_array_indices=0, num_ref_images=1, camera_channels=0, num_cores=1,
-                                  reference_property='sum_col_val',
+                                  reference_property=reference_property,
                                   average_images=False, solver=solver, weighting_preference=-6e-4,
                                   weighting_curvature=1e-7,
                                   num_iterations=2000, lambda_reg=1e-3)
@@ -331,15 +334,36 @@ def create_cc_matrix_file(self):
         file.write("2,3,4\n1,2,7\n3,4,5")
         file.close()
 
-def load_extinction_coefficients_computed(solver, channel, led_array):
-    filename = f'extinction_coefficients_{solver}_channel_{channel}_sum_col_val_led_array_{led_array}.csv'
+    @keyword
+    def strip_bgsub_column_from_led_position_files(self, channel=0):
+        """Remove the bgsub_sum_col_value column from the step 3 output files to simulate
+        data extracted with a ledsa version prior to the local background subtraction."""
+        pattern = os.path.join('analysis', f'channel{channel}', '*_led_positions.csv')
+        for path in glob.glob(pattern):
+            with open(path) as file:
+                lines = file.readlines()
+            with open(path, 'w') as file:
+                for line in lines:
+                    if not line.strip() or line.lstrip().startswith('#'):
+                        file.write(line)
+                    else:
+                        file.write(','.join(line.strip().split(',')[:5]) + '\n')
+
+    @keyword
+    def check_extinction_coefficient_files_are_identical(self, file_a, file_b):
+        data_a = np.loadtxt(file_a, delimiter=',')
+        data_b = np.loadtxt(file_b, delimiter=',')
+        return bool(np.array_equal(data_a, data_b))
+
+def load_extinction_coefficients_computed(solver, channel, led_array, reference_property='sum_col_val'):
+    filename = f'extinction_coefficients_{solver}_channel_{channel}_{reference_property}_led_array_{led_array}.csv'
     data = np.loadtxt(
         os.path.join('analysis', 'extinction_coefficients', solver, filename),delimiter=',')
     time = data[:, 0]
     extinction_coefficients_computed = data[:, 1:]
     return time, extinction_coefficients_computed
 
-def create_test_image(image_id, experiment, img_dir='test_data'):
+def create_test_image(image_id, experiment, img_dir='test_data', background=0):
     """ Creates three test images with black and gray pixels representing 3 leds and sets the exif data needed
     The first image has 100% transmission on all LEDs, the second image has 50% transmission on all LEDs,
     the third has 50%, 70% and 80% transmission on the top, middle and bottom LEDs.
@@ -353,7 +377,7 @@ def create_test_image(image_id, experiment, img_dir='test_data'):
     transmissions = experiment.calc_all_led_transmissions()
 
     # Reverse transmissions because images are created from top down
-    img_array = create_img_array(num_of_leds, list(reversed(transmissions)))
+    img_array = create_img_array(num_of_leds, list(reversed(transmissions)), background=background)
     img = Image.fromarray(img_array, 'RGB')
 
     # Save image without EXIF data
@@ -369,10 +393,12 @@ def create_test_image(image_id, experiment, img_dir='test_data'):
     img2.writeMetadata()
 
 
-def create_img_array(num_of_leds, transmissions):
+def create_img_array(num_of_leds, transmissions, background=0):
     img = np.zeros((num_of_leds * 50 + 50, 50, 3), np.uint8)
     for led_id in range(num_of_leds):
         add_led(img, (1 + led_id) * 50, 25, transmissions[led_id])
+    if background:
+        img = np.clip(img.astype(np.int16) + int(background), 0, 255).astype(np.uint8)
     return img
 
 

ledsa/tests/UnitTests/test_background_subtraction.py

@@ -0,0 +1,118 @@
+import numpy as np
+import pytest
+
+from ledsa.data_extraction.step_3_functions import _estimate_local_background, _generate_led_analysis_data
+
+SEARCH_AREA_RADIUS = 10
+LED_CENTER = (30, 30)
+LED_AMPLITUDE = 200.0
+LED_SIGMA = 2.0
+
+
+def make_synthetic_image(background=0.0, bayer_mask=False):
+    """Create a 60x60 image with a Gaussian LED spot at LED_CENTER plus a constant background.
+
+    With bayer_mask=True, three out of four pixels are set to 0 like in the Bayer
+    array of a raw image where all pixels of foreign color channels are masked with 0.
+    """
+    x, y = np.meshgrid(np.arange(60), np.arange(60), indexing='ij')
+    dist_sq = (x - LED_CENTER[0]) ** 2 + (y - LED_CENTER[1]) ** 2
+    led = LED_AMPLITUDE * np.exp(-dist_sq / (2 * LED_SIGMA ** 2))
+    img = led + background
+    if bayer_mask:
+        mask = (x % 2 == 0) & (y % 2 == 0)
+        img = np.where(mask, img, 0.0)
+    return img
+
+
+def pure_led_integral(img_with_zero_background):
+    radius = SEARCH_AREA_RADIUS
+    return np.sum(img_with_zero_background[LED_CENTER[0] - radius:LED_CENTER[0] + radius,
+                                           LED_CENTER[1] - radius:LED_CENTER[1] + radius])
+
+
+def extract_led_data(img):
+    search_areas = np.array([[0, LED_CENTER[0], LED_CENTER[1]]])
+    return _generate_led_analysis_data(None, 0, img, False, 0, 'test_img', 0, search_areas,
+                                       SEARCH_AREA_RADIUS, fit_leds=False)
+
+
+class TestEstimateLocalBackground:
+    def test_constant_image_returns_background(self):
+        img = np.full((20, 20), 7.5)
+        assert _estimate_local_background(img) == pytest.approx(7.5)
+
+    def test_masked_zero_pixels_are_ignored(self):
+        img = np.full((20, 20), 7.5)
+        img[::2, :] = 0.0
+        img[:, ::2] = 0.0
+        assert _estimate_local_background(img) == pytest.approx(7.5)
+
+    def test_all_zero_border_returns_zero(self):
+        img = np.zeros((20, 20))
+        img[10, 10] = 100.0
+        assert _estimate_local_background(img) == 0.0
+
+    def test_too_small_search_area_returns_zero(self):
+        img = np.full((4, 4), 7.5)
+        assert _estimate_local_background(img) == 0.0
+
+
+class TestBackgroundSubtractedSum:
+    def test_bgsub_sum_is_independent_of_background(self):
+        led_integral = pure_led_integral(make_synthetic_image(background=0.0))
+        for background in [0.0, 5.0, 20.0, 50.0]:
+            led_data = extract_led_data(make_synthetic_image(background=background))
+            assert led_data.bgsub_sum_color_value == pytest.approx(led_integral, rel=1e-3), \
+                f'bgsub_sum_color_value deviates from the pure LED integral for background {background}'
+
+    def test_sum_col_value_scales_with_background(self):
+        led_integral = pure_led_integral(make_synthetic_image(background=0.0))
+        num_pixels = (2 * SEARCH_AREA_RADIUS) ** 2
+        for background in [5.0, 20.0]:
+            led_data = extract_led_data(make_synthetic_image(background=background))
+            assert led_data.sum_color_value == pytest.approx(led_integral + background * num_pixels, rel=1e-6)
+
+    def test_bgsub_sum_with_bayer_masked_pixels(self):
+        led_integral = pure_led_integral(make_synthetic_image(background=0.0, bayer_mask=True))
+        for background in [5.0, 20.0]:
+            led_data = extract_led_data(make_synthetic_image(background=background, bayer_mask=True))
+            assert led_data.bgsub_sum_color_value == pytest.approx(led_integral, rel=1e-3)
+
+    def test_sum_col_value_is_unchanged_by_new_quantity(self):
+        img = make_synthetic_image(background=10.0)
+        led_data = extract_led_data(img)
+        radius = SEARCH_AREA_RADIUS
+        expected = np.sum(img[LED_CENTER[0] - radius:LED_CENTER[0] + radius,
+                              LED_CENTER[1] - radius:LED_CENTER[1] + radius])
+        assert led_data.sum_color_value == pytest.approx(expected)
+
+    def test_csv_serialization_contains_bgsub_value(self):
+        led_data = extract_led_data(make_synthetic_image(background=10.0))
+        main_data_fields = led_data.get_main_data_string().split(',')
+        assert len(main_data_fields) == 6
+        assert float(main_data_fields[5]) == pytest.approx(led_data.bgsub_sum_color_value, rel=1e-4)
+
+
+class TestColumnNames:
+    def _write_led_positions_csv(self, tmp_path, num_columns):
+        channel_dir = tmp_path / 'analysis' / 'channel0'
+        channel_dir.mkdir(parents=True)
+        row = ','.join(str(float(i)) for i in range(num_columns))
+        (channel_dir / '1_led_positions.csv').write_text(f'# header\n{row}\n{row}\n')
+
+    def test_new_csv_with_bgsub_column(self, tmp_path, monkeypatch):
+        from ledsa.core.file_handling import _get_column_names
+        self._write_led_positions_csv(tmp_path, 6)
+        monkeypatch.chdir(tmp_path)
+        columns = _get_column_names(0)
+        assert columns == ["img_id", "led_id", "led_array_id", "sum_col_val", "mean_col_val",
+                           "max_col_val", "bgsub_sum_col_val", "width", "height"]
+
+    def test_old_csv_without_bgsub_column(self, tmp_path, monkeypatch):
+        from ledsa.core.file_handling import _get_column_names
+        self._write_led_positions_csv(tmp_path, 5)
+        monkeypatch.chdir(tmp_path)
+        columns = _get_column_names(0)
+        assert columns == ["img_id", "led_id", "led_array_id", "sum_col_val", "mean_col_val",
+                           "max_col_val", "width", "height"]