diff --git a/ledsa/analysis/ConfigDataAnalysis.py b/ledsa/analysis/ConfigDataAnalysis.py index d8f48b2..ac1d12e 100644 --- a/ledsa/analysis/ConfigDataAnalysis.py +++ b/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 diff --git a/ledsa/core/file_handling.py b/ledsa/core/file_handling.py index 9f7aad9..0fa40c2 100644 --- a/ledsa/core/file_handling.py +++ b/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"]) diff --git a/ledsa/data_extraction/LEDAnalysisData.py b/ledsa/data_extraction/LEDAnalysisData.py index 028706e..d2d43d3 100644 --- a/ledsa/data_extraction/LEDAnalysisData.py +++ b/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: diff --git a/ledsa/data_extraction/step_3_functions.py b/ledsa/data_extraction/step_3_functions.py index e874f9b..d8c2460 100644 --- a/ledsa/data_extraction/step_3_functions.py +++ b/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" diff --git a/ledsa/tests/AcceptanceTests/07_test_background_subtraction.robot b/ledsa/tests/AcceptanceTests/07_test_background_subtraction.robot new file mode 100644 index 0000000..762d709 --- /dev/null +++ b/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 diff --git a/ledsa/tests/AcceptanceTests/LedsaATestLibrary.py b/ledsa/tests/AcceptanceTests/LedsaATestLibrary.py index d616b12..6ca7375 100644 --- a/ledsa/tests/AcceptanceTests/LedsaATestLibrary.py +++ b/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 diff --git a/ledsa/tests/UnitTests/__init__.py b/ledsa/tests/UnitTests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/ledsa/tests/UnitTests/test_background_subtraction.py b/ledsa/tests/UnitTests/test_background_subtraction.py new file mode 100644 index 0000000..011b124 --- /dev/null +++ b/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"]