Take the black level of RAW images from the optical black

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)
 

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

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:

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)

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]

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