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/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_file_handling.py b/ledsa/tests/UnitTests/test_file_handling.py index 63d5548..55f015b 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)) + 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 + + 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] 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