From d644eee3f482ec5c2c204fd27a03caaee6d37261 Mon Sep 17 00:00:00 2001 From: Kristian Date: Tue, 27 Oct 2026 23:17:18 +0100 Subject: [PATCH 1/3] take the black level of RAW images from the optical black The black level of the RAW metadata is a nominal value of the camera model. The optical black pixels are masked pixels at the margin of the sensor, read out together with the image itself, so their mean is the black level of this image and follows the drift with the sensor temperature. The two differ by a fraction of a count, and the integration over the search area of an LED turns that difference into a constant offset of its intensity, which stays while the LED dims in smoke and so saturates the extinction coefficients. _black_level resolves this for the only place that subtracts a black level, so it reaches sum_col_val, bgsub_sum_col_val, sidebg_sum_col_val and bgsub2_sum_col_val alike. The metadata value remains the fallback for sensors whose optical black margin is narrower than MIN_OPTICAL_BLACK_WIDTH. JPG and PNG have no black level and are untouched. The values of the four columns above therefore change for RAW images, which the column layout alone does not show, so the step 3 header now records the LEDSA version as the headers of the analysis results already do. The change only reaches existing experiments when step 3 is run again for them. The fake RAW image of the orientation tests gains the margin sizes and so covers the fallback. Unit tests pass, 185, with test_black_level.py for the estimate and the fallback. Acceptance tests not run, and no RAW image is covered by the unit tests. Co-Authored-By: Claude Opus 5 --- ledsa/core/image_reading.py | 49 ++++++++++++- ledsa/data_extraction/step_3_functions.py | 4 +- ledsa/tests/UnitTests/test_black_level.py | 70 +++++++++++++++++++ .../tests/UnitTests/test_image_orientation.py | 4 +- 4 files changed, 124 insertions(+), 3 deletions(-) create mode 100644 ledsa/tests/UnitTests/test_black_level.py diff --git a/ledsa/core/image_reading.py b/ledsa/core/image_reading.py index 714b7e9..1cd0af8 100644 --- a/ledsa/core/image_reading.py +++ b/ledsa/core/image_reading.py @@ -244,13 +244,60 @@ def _read_channel_data_from_raw_file(filename: str, channel: int, orientation: O return channel_array +# Minimum width of the optical black margin in pixels to estimate the black level from it +MIN_OPTICAL_BLACK_WIDTH = 8 + + +def _optical_black_mean(raw) -> Optional[np.ndarray]: + """ + Mean of the optical black pixels of a RAW image for every color of the Bayer pattern (index as in + raw_colors). The optical black pixels are the masked pixels at the left margin of the sensor, which is + read out with the image itself. + + :param raw: The opened RAW image (rawpy). + :return: Mean per Bayer color, None if the sensor has no optical black margin of at least + MIN_OPTICAL_BLACK_WIDTH pixels. + :rtype: np.ndarray or None + """ + sizes = raw.sizes + if sizes.left_margin < MIN_OPTICAL_BLACK_WIDTH: + return None + rows = np.s_[sizes.top_margin:sizes.top_margin + sizes.height] + margin = raw.raw_image[rows, :sizes.left_margin].astype(np.float64) + colors = raw.raw_colors[rows, :sizes.left_margin] + num_colors = int(raw.raw_colors_visible.max()) + 1 + return np.array([margin[colors == color].mean() for color in range(num_colors)]) + + +def _black_level(raw, channel: int) -> float: + """ + Black level of a color channel, taken from the optical black pixels of the image itself. + + The black level of the RAW metadata is a nominal value of the camera model and is used only as a + fallback, for sensors without an optical black margin. The two differ by a fraction of a count, and the + integration over the search area of an LED turns that difference into a constant offset of its + intensity, which does not vanish when the LED dims in smoke. The mean of the optical black pixels also + follows the drift of the black level with the sensor temperature. + + :param raw: The opened RAW image (rawpy). + :param channel: The color channel to process. + :type channel: int + :return: Black level of the channel. + :rtype: float + """ + optical_black_mean = _optical_black_mean(raw) + if optical_black_mean is None: + return float(raw.black_level_per_channel[channel]) + return float(optical_black_mean[channel]) + + def _read_img_array_from_raw_file(filename: str, channel: int, orientation: Optional[int] = None) -> Tuple[np.ndarray, np.ndarray]: # TODO: channel is only relevant for black level, consider individually! with rawpy.imread(filename) as raw: data = raw.raw_image_visible.copy() filter_array = raw.raw_colors_visible - black_level = raw.black_level_per_channel[channel] + black_level = _black_level(raw, channel) white_level = raw.white_level raw_orientation = RAW_FLIP_ORIENTATIONS.get(raw.sizes.flip, UPRIGHT_ORIENTATION) img_array = data.astype(np.int16) - black_level diff --git a/ledsa/data_extraction/step_3_functions.py b/ledsa/data_extraction/step_3_functions.py index 2cd0af6..f7748bd 100644 --- a/ledsa/data_extraction/step_3_functions.py +++ b/ledsa/data_extraction/step_3_functions.py @@ -1,6 +1,7 @@ import os import re import time +from importlib.metadata import version from typing import List, Optional, Tuple import numpy as np @@ -454,7 +455,8 @@ 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 += f"time[s] = {img_infos[int(img_id) - 1][3]}, " + out_str += f"LEDSA {version('ledsa')}\n" out_str += "# id,line,sum_col_value,average_col_value,max_col_value,bgsub_sum_col_value" out_str += ",sidebg_sum_col_value,sidebg_pedestal,bgsub2_sum_col_value" if fit_leds: diff --git a/ledsa/tests/UnitTests/test_black_level.py b/ledsa/tests/UnitTests/test_black_level.py new file mode 100644 index 0000000..1e49764 --- /dev/null +++ b/ledsa/tests/UnitTests/test_black_level.py @@ -0,0 +1,70 @@ +from types import SimpleNamespace + +import numpy as np +import pytest + +from ledsa.core.image_reading import MIN_OPTICAL_BLACK_WIDTH, _black_level, _optical_black_mean + +# RGGB pattern as in rawpy raw_colors: 0 red, 1 and 3 green, 2 blue +BAYER_TILE = np.array([[0, 1], [3, 2]]) +# black level of every Bayer color in the optical black margin, far from the nominal value below +OPTICAL_BLACK_LEVELS = [100.0, 101.0, 102.0, 103.0] +METADATA_BLACK_LEVELS = [2047, 2047, 2047, 2047] + + +def make_raw(left_margin: int) -> SimpleNamespace: + """ + The part of a rawpy RawPy object that the black level estimate reads: 4 x 4 visible pixels behind a top + and a left margin, with a known black level per Bayer color in the margin. + """ + top_margin, height, width = 2, 4, 4 + shape = (top_margin + height, left_margin + width) + colors = np.tile(BAYER_TILE, (shape[0] // 2 + 1, shape[1] // 2 + 1))[:shape[0], :shape[1]] + raw_image = np.zeros(shape, dtype=np.float64) + for color, level in enumerate(OPTICAL_BLACK_LEVELS): + raw_image[colors == color] = level + return SimpleNamespace( + sizes=SimpleNamespace(left_margin=left_margin, top_margin=top_margin, height=height, width=width), + raw_image=raw_image, + raw_colors=colors, + raw_colors_visible=colors[top_margin:, left_margin:], + black_level_per_channel=METADATA_BLACK_LEVELS, + white_level=16383, + ) + + +class TestOpticalBlackMean: + def test_mean_per_bayer_color(self): + assert _optical_black_mean(make_raw(MIN_OPTICAL_BLACK_WIDTH)) == pytest.approx(OPTICAL_BLACK_LEVELS) + + def test_only_the_margin_enters_the_mean(self): + # the visible pixels of a bright image must not raise the black level + raw = make_raw(MIN_OPTICAL_BLACK_WIDTH) + raw.raw_image[:, MIN_OPTICAL_BLACK_WIDTH:] = 9000.0 + assert _optical_black_mean(raw) == pytest.approx(OPTICAL_BLACK_LEVELS) + + def test_margin_narrower_than_the_minimum_is_not_used(self): + assert _optical_black_mean(make_raw(MIN_OPTICAL_BLACK_WIDTH - 2)) is None + + def test_sensor_without_a_margin(self): + assert _optical_black_mean(make_raw(0)) is None + + +class TestBlackLevel: + @pytest.mark.parametrize('channel', [0, 1, 2]) + def test_optical_black_is_preferred_over_the_metadata(self, channel): + # the metadata hold a nominal value of the camera model, the optical black pixels are read out with + # the image itself + assert _black_level(make_raw(MIN_OPTICAL_BLACK_WIDTH), channel) == \ + pytest.approx(OPTICAL_BLACK_LEVELS[channel]) + + @pytest.mark.parametrize('channel', [0, 1, 2]) + def test_metadata_is_the_fallback_without_an_optical_black_margin(self, channel): + assert _black_level(make_raw(MIN_OPTICAL_BLACK_WIDTH - 1), channel) == \ + pytest.approx(float(METADATA_BLACK_LEVELS[channel])) + + def test_the_black_level_is_a_float_in_both_cases(self): + # it is subtracted from the image and enters the scale factor, where an int would round the + # fraction of a count away that the estimate is about + assert isinstance(_black_level(make_raw(MIN_OPTICAL_BLACK_WIDTH), 0), float) + assert isinstance(_black_level(make_raw(0), 0), float) diff --git a/ledsa/tests/UnitTests/test_image_orientation.py b/ledsa/tests/UnitTests/test_image_orientation.py index 3ddc859..1a69f05 100644 --- a/ledsa/tests/UnitTests/test_image_orientation.py +++ b/ledsa/tests/UnitTests/test_image_orientation.py @@ -163,7 +163,9 @@ class FakeRaw: def __init__(self, raw_image_visible, filter_array, flip): self.raw_image_visible = raw_image_visible self.raw_colors_visible = filter_array - self.sizes = type('Sizes', (), {'flip': flip}) + # a sensor without an optical black margin, so the black level comes from the metadata + self.sizes = type('Sizes', (), {'flip': flip, 'left_margin': 0, 'top_margin': 0, + 'height': raw_image_visible.shape[0]}) self.black_level_per_channel = [0, 0, 0, 0] self.white_level = 255 From 9cf9d2382a3894c85a2bdd8f7c39caca86970e82 Mon Sep 17 00:00:00 2001 From: Kristian Date: Tue, 27 Oct 2026 22:40:35 +0100 Subject: [PATCH 2/3] rewrite the binary data instead of appending to the stale file create_binary_data wrote all_parameters.h5 with append=True and nothing removed the file of an earlier run, so every LED ended up in the table once per run and the analysis saw each image as often as the step was repeated. The binary data are derived from the step 3 result files, so they are now rewritten from scratch and the replacement is reported. The column names are taken from the first result file of a channel and used for all of them, so result files of different width cannot be loaded together. That happens when step 3 is run again for part of the images after a LEDSA version that writes further columns was installed. Until now pandas stopped with "Shape of passed values is (n, 12), indices imply (n, 10)", which says nothing about the cause. The width of every file is now compared with the one of the first file and the error names both files, both counts and the re-run of the data extraction. Unit tests pass, 133, with six new ones in test_file_handling.py. They write a coordinate file, so they take the same path through _append_coordinates as a real analysis. Co-Authored-By: Claude Opus 5 --- ledsa/core/file_handling.py | 22 ++++- ledsa/tests/UnitTests/test_file_handling.py | 104 ++++++++++++++++++++ 2 files changed, 125 insertions(+), 1 deletion(-) diff --git a/ledsa/core/file_handling.py b/ledsa/core/file_handling.py index 0a3bb0b..272e68c 100644 --- a/ledsa/core/file_handling.py +++ b/ledsa/core/file_handling.py @@ -233,7 +233,11 @@ def create_binary_data(channel: int, path: str = '.') -> None: try: os.chdir(path) config = ConfigData() + first_file_path = _find_first_result_file(channel) columns = _get_column_names(channel) + # the names carry the img_id that is prepended on load and the coordinates that are appended below, + # the result file itself has neither + num_columns_of_first_file = len(columns) - 3 fit_params_list = [] @@ -261,6 +265,16 @@ def create_binary_data(channel: int, path: str = '.') -> None: skipped_img_ids.append(image_id) continue + # The column names are taken from the first result file and used for all of them, so files of + # different width cannot be loaded together. They occur when step 3 was run again for part of + # the images after a LEDSA version that writes further columns was installed. + if parameters.shape[-1] != num_columns_of_first_file: + raise ValueError(f'{in_file_path} has {parameters.shape[-1]} columns, while ' + f'{first_file_path} has {num_columns_of_first_file} columns. The result ' + f'files of channel {channel} were written with different column layouts, ' + f'so the column names of the first file do not fit the others. Run the ' + f'data extraction (-s3 or -s3_fast) again for every image of the channel.') + parameters = parameters[parameters[:, 0].argsort()] # sort for led_id parameters = _append_coordinates(parameters) fit_params_fragment = _param_array_to_dataframe(parameters, image_id, columns) @@ -282,7 +296,13 @@ def create_binary_data(channel: int, path: str = '.') -> None: fit_params['max_col_val'] = fit_params['max_col_val'].astype(int) fit_params['sum_col_val'] = fit_params['sum_col_val'].astype(int) out_file_path = os.path.join('analysis', f'channel{channel}', 'all_parameters.h5') - fit_params.to_hdf(out_file_path, key='channel_values', format='table', append=True) + if os.path.exists(out_file_path): + # The binary data are derived from the step 3 files, so they are rewritten from scratch. + # Appending would keep every LED of the earlier run next to the new one, and after a change of + # the columns it fails in the middle of the write. + print(f'{out_file_path} of an earlier run is replaced.') + os.remove(out_file_path) + fit_params.to_hdf(out_file_path, key='channel_values', format='table') finally: os.chdir(prev_dir) diff --git a/ledsa/tests/UnitTests/test_file_handling.py b/ledsa/tests/UnitTests/test_file_handling.py index 63d5548..67d12e5 100644 --- a/ledsa/tests/UnitTests/test_file_handling.py +++ b/ledsa/tests/UnitTests/test_file_handling.py @@ -54,3 +54,107 @@ def test_names_are_read_from_the_first_existing_file(self, tmp_path, monkeypatch monkeypatch.chdir(tmp_path) assert _get_column_names(0) == ["img_id", "led_id", "led_array_id", "sum_col_val", "mean_col_val", "max_col_val", "bgsub_sum_col_val", "width", "height"] + + + +class TestCreateBinaryData: + # led_id, led_array_id, sum_col_val, mean_col_val, max_col_val, bgsub_sum_col_val + NUM_COLUMNS = 6 + # the same file written by v0.11.0, which adds sidebg_sum_col_val, sidebg_pedestal and bgsub2_sum_col_val + NUM_COLUMNS_V011 = 9 + NUM_LEDS = 2 + + @pytest.fixture(autouse=True) + def clean_coordinate_cache(self): + # _append_coordinates caches the coordinate file in a function attribute for the whole process + from ledsa.core.file_handling import _append_coordinates + _append_coordinates.__dict__.pop('coord', None) + yield + _append_coordinates.__dict__.pop('coord', None) + + @pytest.fixture + def analysis_dir(self, tmp_path, monkeypatch): + # create_binary_data only reads the image range out of the config, so a mapping is enough and no + # config.ini has to be written + from ledsa.core import file_handling + monkeypatch.setattr(file_handling, 'ConfigData', + lambda: {'DEFAULT': {'num_img_overflow': 'None'}, + 'analyse_photo': {'first_img_analysis_id': '1', + 'last_img_analysis_id': '2', + 'num_skip_imgs': '0'}}) + channel_dir = tmp_path / 'analysis' / 'channel0' + channel_dir.mkdir(parents=True) + coordinates = '\n'.join(f'{float(led_id)},{float(led_id)},{float(led_id) + 0.5}' + for led_id in range(self.NUM_LEDS)) + (tmp_path / 'analysis' / 'led_search_areas_with_coordinates.csv').write_text(f'{coordinates}\n', + encoding='utf-8') + monkeypatch.chdir(tmp_path) + return tmp_path + + def _write_result_file(self, tmp_path, img_id, num_columns): + rows = '\n'.join(','.join(str(float(led_id * 10 + column)) for column in range(num_columns)) + for led_id in range(self.NUM_LEDS)) + path = tmp_path / 'analysis' / 'channel0' / f'{img_id}_led_positions.csv' + path.write_text(f'# header\n{rows}\n', encoding='utf-8') + + def _load(self, tmp_path): + import pandas as pd + return pd.read_hdf(tmp_path / 'analysis' / 'channel0' / 'all_parameters.h5', key='channel_values') + + def test_every_led_of_every_image_is_loaded(self, analysis_dir): + from ledsa.core.file_handling import create_binary_data + for img_id in (1, 2): + self._write_result_file(analysis_dir, img_id, self.NUM_COLUMNS) + create_binary_data(0) + table = self._load(analysis_dir) + assert len(table) == 2 * self.NUM_LEDS + assert sorted(table['img_id'].unique()) == [1, 2] + + def test_a_binary_file_of_an_earlier_run_is_replaced(self, analysis_dir, capsys): + # appending kept every LED of the earlier run next to the new one, so the rows piled up with + # every run and the analysis saw each image twice + from ledsa.core.file_handling import create_binary_data + for img_id in (1, 2): + self._write_result_file(analysis_dir, img_id, self.NUM_COLUMNS) + create_binary_data(0) + rows_of_the_first_run = len(self._load(analysis_dir)) + create_binary_data(0) + assert len(self._load(analysis_dir)) == rows_of_the_first_run + assert 'is replaced' in capsys.readouterr().out + + def test_a_replaced_binary_file_holds_the_new_values(self, analysis_dir): + from ledsa.core.file_handling import create_binary_data + for img_id in (1, 2): + self._write_result_file(analysis_dir, img_id, self.NUM_COLUMNS) + create_binary_data(0) + for img_id in (1, 2): + self._write_result_file(analysis_dir, img_id, self.NUM_COLUMNS) + rows = self._load(analysis_dir) + create_binary_data(0) + assert self._load(analysis_dir)['sum_col_val'].tolist() == rows['sum_col_val'].tolist() + + def test_differing_column_counts_name_the_re_run(self, analysis_dir): + # step 3 was run again for part of the images with a version that writes further columns + from ledsa.core.file_handling import create_binary_data + self._write_result_file(analysis_dir, 1, self.NUM_COLUMNS) + self._write_result_file(analysis_dir, 2, self.NUM_COLUMNS_V011) + with pytest.raises(ValueError, match=r'2_led_positions\.csv has 9 columns, while .*' + r'1_led_positions\.csv has 6 columns.*' + r'Run the data extraction'): + create_binary_data(0) + + def test_differing_column_counts_are_found_against_the_wider_first_file(self, analysis_dir): + # the names come from the lowest image ID, so the wider file may be the one that is read first + from ledsa.core.file_handling import create_binary_data + self._write_result_file(analysis_dir, 1, self.NUM_COLUMNS_V011) + self._write_result_file(analysis_dir, 2, self.NUM_COLUMNS) + with pytest.raises(ValueError, match=r'2_led_positions\.csv has 6 columns, while .*' + r'1_led_positions\.csv has 9 columns'): + create_binary_data(0) + + def test_an_image_without_results_is_still_skipped(self, analysis_dir): + # image 1 could not be read, so it has no result file and the names come from image 2 + from ledsa.core.file_handling import create_binary_data + self._write_result_file(analysis_dir, 2, self.NUM_COLUMNS) + create_binary_data(0) + assert sorted(self._load(analysis_dir)['img_id'].unique()) == [2] From b37330796bc6df48184e9eb4027ab0dd8cadb857 Mon Sep 17 00:00:00 2001 From: Kristian Date: Tue, 27 Oct 2026 23:22:52 +0100 Subject: [PATCH 3/3] name the encoding on the line of the call The CI check for text I/O without an encoding is line based, so the encoding of a wrapped call reads as missing to it although Python sees it. Co-Authored-By: Claude Opus 5 --- ledsa/tests/UnitTests/test_file_handling.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ledsa/tests/UnitTests/test_file_handling.py b/ledsa/tests/UnitTests/test_file_handling.py index 67d12e5..55f015b 100644 --- a/ledsa/tests/UnitTests/test_file_handling.py +++ b/ledsa/tests/UnitTests/test_file_handling.py @@ -86,8 +86,8 @@ def analysis_dir(self, tmp_path, monkeypatch): channel_dir.mkdir(parents=True) coordinates = '\n'.join(f'{float(led_id)},{float(led_id)},{float(led_id) + 0.5}' for led_id in range(self.NUM_LEDS)) - (tmp_path / 'analysis' / 'led_search_areas_with_coordinates.csv').write_text(f'{coordinates}\n', - encoding='utf-8') + coordinates_path = tmp_path / 'analysis' / 'led_search_areas_with_coordinates.csv' + coordinates_path.write_text(f'{coordinates}\n', encoding='utf-8') monkeypatch.chdir(tmp_path) return tmp_path