From 3fff82611d916172452697a55f96a350d50c18d7 Mon Sep 17 00:00:00 2001 From: Kristian Date: Wed, 28 Oct 2026 22:31:38 +0100 Subject: [PATCH 1/3] read the RAW data without the clipping at zero as well The channel data are clipped at zero after the black level subtraction, which turns the negative half of the read noise into a positive offset of the integrated intensity of every LED. The offset stays while the LED dims in smoke, so it saturates the extinction coefficients in dense smoke, where it is 0.3 to 0.7 % of I0 in the Warsaw tests. _read_img_array_and_optical_black_from_raw_file returns the same data as before plus a version without the clipping, in which every Bayer color keeps the black level of its own optical black pixels, and the read noise of every single pixel from the spread of those pixels. The file is read once for all of them, and the orientation reaches all four arrays. _read_raw_image, _clipped_data and _unclipped_data split what was one function, so the two readers share the way the data are scaled instead of repeating it. _optical_black_mean becomes _optical_black_statistics, since the noise needs the spread of the margin next to its mean. Unit tests pass, 315. Co-Authored-By: Claude Opus 5 --- ledsa/core/image_reading.py | 110 ++++++++++++++++---- ledsa/tests/UnitTests/test_black_level.py | 121 +++++++++++++++++++--- 2 files changed, 191 insertions(+), 40 deletions(-) diff --git a/ledsa/core/image_reading.py b/ledsa/core/image_reading.py index 1cd0af8..e006d6d 100644 --- a/ledsa/core/image_reading.py +++ b/ledsa/core/image_reading.py @@ -1,6 +1,6 @@ import os from datetime import datetime, timedelta -from typing import Callable, Optional, Tuple +from typing import Callable, NamedTuple, Optional, Tuple import exiv2 import numpy as np @@ -248,25 +248,28 @@ def _read_channel_data_from_raw_file(filename: str, channel: int, orientation: O MIN_OPTICAL_BLACK_WIDTH = 8 -def _optical_black_mean(raw) -> Optional[np.ndarray]: +def _optical_black_statistics(raw) -> Tuple[Optional[np.ndarray], 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. + Mean and standard deviation 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, so the mean is the black level of this image and + the standard deviation is the read noise of a single pixel. :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 + :return: Mean and standard deviation per Bayer color, (None, None) if the sensor has no optical black + margin of at least MIN_OPTICAL_BLACK_WIDTH pixels. + :rtype: tuple[np.ndarray or None, np.ndarray or None] """ sizes = raw.sizes if sizes.left_margin < MIN_OPTICAL_BLACK_WIDTH: - return None + return None, 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)]) + per_color = [margin[colors == color] for color in range(num_colors)] + return (np.array([values.mean() for values in per_color]), + np.array([values.std() for values in per_color])) def _black_level(raw, channel: int) -> float: @@ -285,28 +288,89 @@ def _black_level(raw, channel: int) -> float: :return: Black level of the channel. :rtype: float """ - optical_black_mean = _optical_black_mean(raw) + optical_black_mean, _ = _optical_black_statistics(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]: +class RawImage(NamedTuple): + """What LEDSA reads from a RAW file, in one pass and before the orientation is applied.""" + + data: np.ndarray # the visible sensor values as they are stored + filter_array: np.ndarray # the Bayer color of every pixel + black_level: float # of the selected channel + white_level: float + orientation: int # the one the camera recorded in its maker notes + optical_black_mean: Optional[np.ndarray] # per Bayer color, None without an optical black margin + optical_black_std: Optional[np.ndarray] # per Bayer color, None without an optical black margin + + @property + def scale(self) -> float: + """Factor that maps the range above the black level onto the full range of the sensor.""" + return self.white_level / (self.white_level - self.black_level) + + +def _read_raw_image(filename: str, channel: int) -> RawImage: # 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 = _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 - img_array = (img_array * (white_level / (white_level - black_level))).astype(np.int16) - img_array = np.clip(img_array, 0, white_level) + optical_black_mean, optical_black_std = _optical_black_statistics(raw) + black_level = (float(raw.black_level_per_channel[channel]) if optical_black_mean is None + else float(optical_black_mean[channel])) + return RawImage(raw.raw_image_visible.copy(), raw.raw_colors_visible.copy(), black_level, + raw.white_level, RAW_FLIP_ORIENTATIONS.get(raw.sizes.flip, UPRIGHT_ORIENTATION), + optical_black_mean, optical_black_std) + + +def _clipped_data(raw: RawImage) -> np.ndarray: + """ + The sensor values with the black level of the channel subtracted, scaled to the white level and clipped + at zero, which is how LEDSA has always read them. + """ + img_array = (raw.data.astype(np.int16) - raw.black_level) * raw.scale + return np.clip(img_array.astype(np.int16), 0, raw.white_level) + + +def _unclipped_data(raw: RawImage) -> np.ndarray: + """ + The same values without the clipping at zero and with the black level of every Bayer color taken + individually. The clipping turns the negative half of the read noise into a positive offset of the + integrated intensity of an LED, which stays while the LED dims in smoke; these values have no such + offset and are negative where the sensor read below its black level. + """ + corrected = raw.data.astype(np.float32) - raw.optical_black_mean.astype(np.float32)[raw.filter_array] + return corrected * np.float32(raw.scale) + + +def _read_img_array_from_raw_file(filename: str, channel: int, + orientation: Optional[int] = None) -> Tuple[np.ndarray, np.ndarray]: + raw = _read_raw_image(filename, channel) if orientation is None: - orientation = raw_orientation + orientation = raw.orientation # The Bayer pattern is turned with the pixels, so the filter array keeps naming the color of every pixel - return _apply_orientation(img_array, orientation), _apply_orientation(filter_array, orientation) + return (_apply_orientation(_clipped_data(raw), orientation), + _apply_orientation(raw.filter_array, orientation)) + + +def _read_img_array_and_optical_black_from_raw_file( + filename: str, channel: int, + orientation: Optional[int] = None) -> Tuple[np.ndarray, np.ndarray, Optional[np.ndarray], + Optional[np.ndarray]]: + """ + As _read_img_array_from_raw_file, plus the data without the clipping at zero and the read noise of every + single pixel. The file is read once for all of them. + """ + raw = _read_raw_image(filename, channel) + if orientation is None: + orientation = raw.orientation + img_array = _apply_orientation(_clipped_data(raw), orientation) + filter_array = _apply_orientation(raw.filter_array, orientation) + if raw.optical_black_mean is None: + return img_array, filter_array, None, None + # the noise of a pixel is the one of the optical black pixels of its own Bayer color + noise = (raw.optical_black_std.astype(np.float32) * np.float32(raw.scale))[raw.filter_array] + return (img_array, filter_array, _apply_orientation(_unclipped_data(raw), orientation), + _apply_orientation(noise, orientation)) def _read_img_array_from_img_file(filename: str, orientation: Optional[int] = None) -> np.ndarray: img_array = plt.imread(filename) diff --git a/ledsa/tests/UnitTests/test_black_level.py b/ledsa/tests/UnitTests/test_black_level.py index 1e49764..f5d1b4e 100644 --- a/ledsa/tests/UnitTests/test_black_level.py +++ b/ledsa/tests/UnitTests/test_black_level.py @@ -3,19 +3,26 @@ import numpy as np import pytest -from ledsa.core.image_reading import MIN_OPTICAL_BLACK_WIDTH, _black_level, _optical_black_mean +import ledsa.core.image_reading as image_reading +from ledsa.core.image_reading import (MIN_OPTICAL_BLACK_WIDTH, UPRIGHT_ORIENTATION, _black_level, + _optical_black_statistics, + _read_img_array_and_optical_black_from_raw_file, + _read_img_array_from_raw_file) # 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] +WHITE_LEVEL = 16383 -def make_raw(left_margin: int) -> SimpleNamespace: +def make_raw(left_margin: int, noise: float = 0.0, visible: float = None) -> 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. + The part of a rawpy RawPy object that LEDSA reads: 4 x 4 visible pixels behind a top and a left margin, + with a known black level per Bayer color in the margin. With noise the two halves of the margin columns + are offset by -noise and +noise, so the standard deviation of every color is noise. With visible the + visible pixels are set to that value instead of their black level. """ top_margin, height, width = 2, 4, 4 shape = (top_margin + height, left_margin + width) @@ -23,31 +30,57 @@ def make_raw(left_margin: int) -> SimpleNamespace: 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), + if noise and left_margin: + raw_image[:, :left_margin] += np.where(np.arange(left_margin) < left_margin // 2, -noise, noise) + if visible is not None: + raw_image[top_margin:, left_margin:] = visible + raw = SimpleNamespace( + sizes=SimpleNamespace(left_margin=left_margin, top_margin=top_margin, height=height, width=width, + flip=0), raw_image=raw_image, + raw_image_visible=raw_image[top_margin:, left_margin:], raw_colors=colors, raw_colors_visible=colors[top_margin:, left_margin:], black_level_per_channel=METADATA_BLACK_LEVELS, - white_level=16383, + white_level=WHITE_LEVEL, ) + raw.__enter__ = lambda: raw + raw.__exit__ = lambda *args: False + return raw -class TestOpticalBlackMean: - def test_mean_per_bayer_color(self): - assert _optical_black_mean(make_raw(MIN_OPTICAL_BLACK_WIDTH)) == pytest.approx(OPTICAL_BLACK_LEVELS) +@pytest.fixture +def fake_imread(monkeypatch): + """Let rawpy.imread return a prepared fake RAW image.""" + def install(raw): + class Opened: + def __enter__(self): + return raw - 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 __exit__(self, *args): + return False + + monkeypatch.setattr(image_reading.rawpy, 'imread', lambda filename: Opened()) + return install + + +class TestOpticalBlackStatistics: + def test_mean_and_standard_deviation_per_bayer_color(self): + mean, std = _optical_black_statistics(make_raw(MIN_OPTICAL_BLACK_WIDTH, noise=3.0)) + assert mean == pytest.approx(OPTICAL_BLACK_LEVELS) + assert std == pytest.approx([3.0] * len(OPTICAL_BLACK_LEVELS)) + + def test_only_the_margin_enters_the_statistics(self): + # the visible pixels of a bright image must not raise the black level or the noise + mean, std = _optical_black_statistics(make_raw(MIN_OPTICAL_BLACK_WIDTH, visible=9000.0)) + assert mean == pytest.approx(OPTICAL_BLACK_LEVELS) + assert std == pytest.approx([0.0] * len(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 + assert _optical_black_statistics(make_raw(MIN_OPTICAL_BLACK_WIDTH - 2)) == (None, None) def test_sensor_without_a_margin(self): - assert _optical_black_mean(make_raw(0)) is None + assert _optical_black_statistics(make_raw(0)) == (None, None) class TestBlackLevel: @@ -68,3 +101,57 @@ def test_the_black_level_is_a_float_in_both_cases(self): # 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) + + +class TestUnclippedData: + CHANNEL = 0 + + def test_values_below_the_black_level_stay_negative(self, fake_imread): + # the visible pixels read 2 counts below the black level of their color. Only the pixels of the + # selected channel are compared, because the clipped version subtracts one black level from the + # whole Bayer array (see the TODO on _read_raw_image). + raw = make_raw(MIN_OPTICAL_BLACK_WIDTH) + raw.raw_image[raw.sizes.top_margin:, MIN_OPTICAL_BLACK_WIDTH:] -= 2.0 + fake_imread(raw) + clipped, filter_array, unclipped, _ = \ + _read_img_array_and_optical_black_from_raw_file('img.CR3', self.CHANNEL) + mask = filter_array == self.CHANNEL + scale = WHITE_LEVEL / (WHITE_LEVEL - OPTICAL_BLACK_LEVELS[self.CHANNEL]) + assert (clipped[mask] == 0).all() # the whole signal is gone, and this is the offset it leaves + assert unclipped[mask] == pytest.approx(np.full(mask.sum(), -2.0 * scale), rel=1e-5) + + def test_every_bayer_color_gets_its_own_black_level(self, fake_imread): + # the visible pixels sit exactly at the black level of their own color, so all of them are 0 + fake_imread(make_raw(MIN_OPTICAL_BLACK_WIDTH)) + _, _, unclipped, _ = _read_img_array_and_optical_black_from_raw_file('img.CR3', self.CHANNEL) + assert unclipped == pytest.approx(np.zeros(unclipped.shape)) + + def test_the_noise_of_a_pixel_is_the_one_of_its_color(self, fake_imread): + fake_imread(make_raw(MIN_OPTICAL_BLACK_WIDTH, noise=4.0)) + _, filter_array, _, noise = _read_img_array_and_optical_black_from_raw_file('img.CR3', self.CHANNEL) + scale = WHITE_LEVEL / (WHITE_LEVEL - OPTICAL_BLACK_LEVELS[self.CHANNEL]) + assert noise.shape == filter_array.shape + assert noise == pytest.approx(np.full(noise.shape, 4.0 * scale), rel=1e-5) + + def test_the_clipped_version_is_the_one_every_other_reader_gets(self, fake_imread): + fake_imread(make_raw(MIN_OPTICAL_BLACK_WIDTH, visible=5000.0)) + clipped, filter_array, _, _ = _read_img_array_and_optical_black_from_raw_file('img.CR3', self.CHANNEL) + expected, expected_filter = _read_img_array_from_raw_file('img.CR3', self.CHANNEL) + np.testing.assert_array_equal(clipped, expected) + np.testing.assert_array_equal(filter_array, expected_filter) + + def test_without_an_optical_black_margin_there_is_nothing_to_return(self, fake_imread): + fake_imread(make_raw(MIN_OPTICAL_BLACK_WIDTH - 1, visible=5000.0)) + clipped, filter_array, unclipped, noise = \ + _read_img_array_and_optical_black_from_raw_file('img.CR3', self.CHANNEL) + assert unclipped is None and noise is None + assert clipped.shape == filter_array.shape + + def test_the_orientation_reaches_all_four_arrays(self, fake_imread): + raw = make_raw(MIN_OPTICAL_BLACK_WIDTH, noise=4.0, visible=5000.0) + fake_imread(raw) + upright = _read_img_array_and_optical_black_from_raw_file('img.CR3', self.CHANNEL, + UPRIGHT_ORIENTATION) + turned = _read_img_array_and_optical_black_from_raw_file('img.CR3', self.CHANNEL, 6) + for upright_array, turned_array in zip(upright, turned): + np.testing.assert_array_equal(turned_array, np.rot90(upright_array, -1)) From 35fc09dae5010f32de18dc78311aa64b095ff61f Mon Sep 17 00:00:00 2001 From: Kristian Date: Wed, 28 Oct 2026 22:34:11 +0100 Subject: [PATCH 2/3] write the unclipped intensity and its noise in step 3 Step 3 writes two columns after the side band ones: ob_sum_col_value, the integrated intensity of an LED without the clipping at zero, and ob_noise_sum_col_value, the read noise of that sum. Both are NaN for JPG and PNG images and for sensors without an optical black margin. sum_col_val and the columns before them keep their place and their values, so the prefix of the result files released with v0.11.0 is unchanged. The layouts of -s3_fast and -s3 now have 12 and 26 names, which collide with none of the earlier ones. read_channel_data_and_mask_from_img becomes read_channel_data_mask_and_optical_black_from_img and returns the two further arrays, so step 3 still reads every image once. Both carry 0 outside the color channel, so the sum over a search area and the variance of that sum need the channel mask no further. NOISE_COLUMNS names the noise column of a reference property. Only ob_sum_col_val has one, which is what the weighting of the next commit needs. The RAW test of the channel mask now mocks the one pass over the RAW file instead of one reader, so every reader derived from it stays consistent in the test. Unit tests pass, 320. Co-Authored-By: Claude Opus 5 --- ledsa/core/file_handling.py | 11 +++- ledsa/core/image_reading.py | 28 +++++++--- ledsa/data_extraction/LEDAnalysisData.py | 10 ++++ ledsa/data_extraction/step_3_functions.py | 20 ++++++- ledsa/tests/UnitTests/test_black_level.py | 56 +++++++++++++++++++ .../tests/UnitTests/test_image_orientation.py | 9 ++- .../UnitTests/test_sideband_background.py | 16 ++++-- 7 files changed, 130 insertions(+), 20 deletions(-) diff --git a/ledsa/core/file_handling.py b/ledsa/core/file_handling.py index cb6415a..0beb6ef 100644 --- a/ledsa/core/file_handling.py +++ b/ledsa/core/file_handling.py @@ -517,6 +517,7 @@ def create_binary_data(channel: int, path: str = '.') -> None: BASE_COLUMNS = ["img_id", "led_id", "led_array_id", "sum_col_val", "mean_col_val", "max_col_val"] BGSUB_COLUMNS = ["bgsub_sum_col_val"] SIDEBG_COLUMNS = ["sidebg_sum_col_val", "sidebg_pedestal", "bgsub2_sum_col_val"] +OB_COLUMNS = ["ob_sum_col_val", "ob_noise_sum_col_val"] LED_CENTER_COLUMNS = ["led_center_x", "led_center_y"] FIT_COLUMNS = ["x", "y", "dx", "dy", "A", "alpha", "wx", "wy", "fit_success", "fit_fun", "fit_nfev", "fit_time"] @@ -528,8 +529,16 @@ def create_binary_data(channel: int, path: str = '.') -> None: len(BASE_COLUMNS + BGSUB_COLUMNS + LED_CENTER_COLUMNS + FIT_COLUMNS): BASE_COLUMNS + BGSUB_COLUMNS + LED_CENTER_COLUMNS + FIT_COLUMNS, # -s3 with fit, before the side band pedestal len(BASE_COLUMNS + BGSUB_COLUMNS + SIDEBG_COLUMNS + LED_CENTER_COLUMNS + FIT_COLUMNS): - BASE_COLUMNS + BGSUB_COLUMNS + SIDEBG_COLUMNS + LED_CENTER_COLUMNS + FIT_COLUMNS, # -s3 with fit + BASE_COLUMNS + BGSUB_COLUMNS + SIDEBG_COLUMNS + LED_CENTER_COLUMNS + FIT_COLUMNS, + # -s3 with fit, before the optical black + len(BASE_COLUMNS + BGSUB_COLUMNS + SIDEBG_COLUMNS + OB_COLUMNS): + BASE_COLUMNS + BGSUB_COLUMNS + SIDEBG_COLUMNS + OB_COLUMNS, # -s3_fast + len(BASE_COLUMNS + BGSUB_COLUMNS + SIDEBG_COLUMNS + OB_COLUMNS + LED_CENTER_COLUMNS + FIT_COLUMNS): + BASE_COLUMNS + BGSUB_COLUMNS + SIDEBG_COLUMNS + OB_COLUMNS + LED_CENTER_COLUMNS + FIT_COLUMNS, + # -s3 with fit } +# noise column of a reference property, used to weight the LEDs in the linear solver +NOISE_COLUMNS = {"ob_sum_col_val": "ob_noise_sum_col_val"} # Fitted files written before the delimiter between the main and the fit part was added: max_col_value and # led_center_x are glued into one field, so the file has one column less and the names are shifted from there on. LEGACY_FIT_COLUMNS = BASE_COLUMNS + LED_CENTER_COLUMNS + FIT_COLUMNS[:-1] diff --git a/ledsa/core/image_reading.py b/ledsa/core/image_reading.py index e006d6d..582af3c 100644 --- a/ledsa/core/image_reading.py +++ b/ledsa/core/image_reading.py @@ -72,16 +72,22 @@ def read_channel_data_from_img(filename: str, channel: int, rotation: str = AUTO return _read_or_raise(filename, reader, channel, orientation) -def read_channel_data_and_mask_from_img(filename: str, channel: int, - rotation: str = AUTO_ROTATION) -> Tuple[np.ndarray, np.ndarray]: +def read_channel_data_mask_and_optical_black_from_img( + filename: str, channel: int, + rotation: str = AUTO_ROTATION) -> Tuple[np.ndarray, np.ndarray, Optional[np.ndarray], + Optional[np.ndarray]]: """ Returns the same 2D channel array as read_channel_data_from_img together with a boolean mask of the - pixels that belong to the selected color channel. + pixels that belong to the selected color channel and, for RAW images with an optical black margin, the + same data without the clipping at zero and the read noise of every pixel. For RAW files the mask is taken from the Bayer filter pattern (raw_colors_visible) and not from the pixel values, because pixels of the selected channel can be 0 as well (noise clipped at 0 after the black level subtraction). For JPG/PNG files every pixel belongs to the channel. + The unclipped data and the noise carry 0 outside the channel, like the channel array itself, so that a + sum over a search area and the noise of that sum need the mask no further. + :param filename: The path of the image file to read. :type filename: str :param channel: The color channel to process. @@ -89,19 +95,25 @@ def read_channel_data_and_mask_from_img(filename: str, channel: int, :param rotation: Degrees the stored array is turned clockwise to be upright, 'auto' to take the rotation from the metadata of the image. :type rotation: str or int - :return: The channel array and the boolean channel mask of the same shape. - :rtype: tuple[np.ndarray, np.ndarray] + :return: The channel array, the boolean channel mask, the unclipped channel array and the read noise of + every pixel, all of the same shape. The last two are None for JPG and PNG images and for sensors + without an optical black margin. + :rtype: tuple[np.ndarray, np.ndarray, np.ndarray or None, np.ndarray or None] :raises CorruptedImageError: If the image cannot be read or has an unsupported file extension. :raises ValueError: If rotation is neither 'auto' nor one of 0, 90, 180 and 270. """ def read_img_file(img_filename, img_channel, img_orientation=None): channel_array = _read_channel_data_from_img_file(img_filename, img_channel, img_orientation) - return channel_array, np.ones(channel_array.shape, dtype=bool) + return channel_array, np.ones(channel_array.shape, dtype=bool), None, None def read_raw_file(raw_filename, raw_channel, raw_orientation=None): - img_array, filter_array = _read_img_array_from_raw_file(raw_filename, raw_channel, raw_orientation) + img_array, filter_array, unclipped, noise = _read_img_array_and_optical_black_from_raw_file( + raw_filename, raw_channel, raw_orientation) channel_mask = get_bayer_channel_mask(filter_array, raw_channel) - return np.where(channel_mask, img_array, 0), channel_mask + if unclipped is None: + return np.where(channel_mask, img_array, 0), channel_mask, None, None + return (np.where(channel_mask, img_array, 0), channel_mask, + np.where(channel_mask, unclipped, 0.0), np.where(channel_mask, noise, 0.0)) orientation = _orientation_from_rotation(rotation) reader = _get_reader(filename, read_img_file, read_raw_file) diff --git a/ledsa/data_extraction/LEDAnalysisData.py b/ledsa/data_extraction/LEDAnalysisData.py index 7c19af8..7ccb331 100644 --- a/ledsa/data_extraction/LEDAnalysisData.py +++ b/ledsa/data_extraction/LEDAnalysisData.py @@ -31,6 +31,13 @@ class LEDAnalysisData: over the channel pixels (including 0) of the border ring of the search area, at least a fraction of the plain sum (background_floor). :vartype bgsub2_sum_color_value: float + :ivar ob_sum_color_value: Integrated color value of the LED over the search area without the clipping at + zero after the black level subtraction, with the black level of every Bayer color taken individually + (RAW images with an optical black margin only, NaN otherwise). + :vartype ob_sum_color_value: float + :ivar ob_noise_sum_color_value: Standard deviation of ob_sum_color_value from the read noise of the + sensor, estimated from the optical black pixels (NaN where ob_sum_color_value is NaN). + :vartype ob_noise_sum_color_value: float :ivar fit_results: Fit results after fitting. :vartype fit_results: OptimizeResult :ivar fit_time: Time taken for fitting. @@ -58,6 +65,8 @@ def __init__(self, led_id, led_array, fit_leds): self.sidebg_sum_color_value = None self.sidebg_pedestal = None self.bgsub2_sum_color_value = None + self.ob_sum_color_value = float('nan') + self.ob_noise_sum_color_value = float('nan') self.fit_leds = fit_leds self.fit_results = None self.fit_time = None @@ -88,6 +97,7 @@ def get_main_data_string(self) -> str: out_str += f',{self.bgsub_sum_color_value:10.4e}' out_str += f',{self.sidebg_sum_color_value:10.4e},{self.sidebg_pedestal:10.4e}' out_str += f',{self.bgsub2_sum_color_value:10.4e}' + out_str += f',{self.ob_sum_color_value:10.4e},{self.ob_noise_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 f7748bd..2e33c8e 100644 --- a/ledsa/data_extraction/step_3_functions.py +++ b/ledsa/data_extraction/step_3_functions.py @@ -10,7 +10,7 @@ from ledsa.core.ConfigData import ConfigData from ledsa.core.file_handling import read_table from ledsa.core.image_handling import get_img_name -from ledsa.core.image_reading import read_channel_data_and_mask_from_img +from ledsa.core.image_reading import read_channel_data_mask_and_optical_black_from_img from ledsa.data_extraction.LEDAnalysisData import LEDAnalysisData from ledsa.data_extraction.model import target_function from ledsa.data_extraction.sideband_background import (DEFAULT_SIDEBAND_DISTANCES, DEFAULT_BACKGROUND_FLOOR, @@ -43,8 +43,8 @@ def generate_analysis_data(img_filename: str, channel: int, search_areas: np.nda :rtype: List[LEDAnalysisData] """ file_path = os.path.join(conf['DEFAULT']['img_directory'], img_filename) - data, channel_mask = read_channel_data_and_mask_from_img(file_path, channel=channel, - rotation=conf.get('DEFAULT', 'img_rotation')) + data, channel_mask, ob_data, ob_noise = read_channel_data_mask_and_optical_black_from_img( + file_path, channel=channel, rotation=conf.get('DEFAULT', 'img_rotation')) search_area_radius = int(conf['find_search_areas']['search_area_radius']) img_analysis_data = [] @@ -64,6 +64,7 @@ def generate_analysis_data(img_filename: str, channel: int, search_areas: np.nda led_analysis_data = _generate_led_analysis_data(conf, channel, data, debug, iled, img_filename, led_array_idx, search_areas, search_area_radius, fit_leds, channel_mask=channel_mask, + ob_data=ob_data, ob_noise=ob_noise, sideband_pixels=sideband_pixels[iled], background_floor=background_floor) img_analysis_data.append(led_analysis_data) @@ -175,6 +176,8 @@ def find_and_save_not_analysed_imgs(channel: int) -> None: def _generate_led_analysis_data(conf: ConfigData, channel: int, data: np.ndarray, debug: bool, iled: int, img_filename: str, led_array_idx: int, search_areas: np.ndarray, search_area_radius: int, fit_leds: bool = True, channel_mask: Optional[np.ndarray] = None, + ob_data: Optional[np.ndarray] = None, + ob_noise: Optional[np.ndarray] = None, sideband_pixels: Optional[Tuple[np.ndarray, np.ndarray]] = None, background_floor: float = DEFAULT_BACKGROUND_FLOOR) -> LEDAnalysisData: """ @@ -200,6 +203,11 @@ def _generate_led_analysis_data(conf: ConfigData, channel: int, data: np.ndarray :type search_area_radius: int :param fit_leds: If True, the LED is fitted to a model function. :type fit_leds: bool + :param ob_data: Channel data without the clipping at zero, None if the image has none. If None, the + optical black values of the LED stay NaN. + :type ob_data: np.ndarray or None + :param ob_noise: Read noise of every pixel of ob_data. + :type ob_noise: np.ndarray or None :param channel_mask: Boolean mask of the pixels belonging to the color channel (Bayer pattern for raw files). None if all pixels belong to the channel. :type channel_mask: np.ndarray or None @@ -238,6 +246,11 @@ def _generate_led_analysis_data(conf: ConfigData, channel: int, data: np.ndarray led_data.bgsub_sum_color_value = float(np.sum(np.clip(data[search_area].astype(np.float64) - background, 0.0, None))) + if ob_data is not None: + # both arrays are 0 outside the channel, so the sum and its variance need no mask + led_data.ob_sum_color_value = float(np.sum(ob_data[search_area], dtype=np.float64)) + led_data.ob_noise_sum_color_value = float(np.sqrt(np.sum(ob_noise[search_area].astype(np.float64) ** 2))) + # the pedestals below are subtracted as N * b with N the number of channel pixels in the box box_mask = None if channel_mask is None else channel_mask[search_area] num_channel_pixels = data[search_area].size if box_mask is None else int(np.count_nonzero(box_mask)) @@ -459,6 +472,7 @@ def _create_header(channel: int, img_id: str, img_filename: str, img_infos: np.n 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" + out_str += ",ob_sum_col_value,ob_noise_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/UnitTests/test_black_level.py b/ledsa/tests/UnitTests/test_black_level.py index f5d1b4e..6db109d 100644 --- a/ledsa/tests/UnitTests/test_black_level.py +++ b/ledsa/tests/UnitTests/test_black_level.py @@ -155,3 +155,59 @@ def test_the_orientation_reaches_all_four_arrays(self, fake_imread): turned = _read_img_array_and_optical_black_from_raw_file('img.CR3', self.CHANNEL, 6) for upright_array, turned_array in zip(upright, turned): np.testing.assert_array_equal(turned_array, np.rot90(upright_array, -1)) + + +class TestOpticalBlackSumsOfAnLed: + """The two columns step 3 writes from the unclipped data.""" + + SEARCH_AREA_RADIUS = 6 + CENTER = (20, 20) + + def _extract(self, data, ob_data=None, ob_noise=None): + from ledsa.data_extraction.step_3_functions import _generate_led_analysis_data + search_areas = np.array([[0, self.CENTER[0], self.CENTER[1]]]) + return _generate_led_analysis_data(None, 0, data, False, 0, 'img', 0, search_areas, + self.SEARCH_AREA_RADIUS, fit_leds=False, + ob_data=ob_data, ob_noise=ob_noise) + + @property + def _box(self): + radius = self.SEARCH_AREA_RADIUS + return np.index_exp[self.CENTER[0] - radius:self.CENTER[0] + radius, + self.CENTER[1] - radius:self.CENTER[1] + radius] + + def test_without_optical_black_the_columns_stay_nan(self): + led_data = self._extract(np.full((40, 40), 5.0)) + assert np.isnan(led_data.ob_sum_color_value) + assert np.isnan(led_data.ob_noise_sum_color_value) + + def test_the_sum_runs_over_the_search_area(self): + data = np.full((40, 40), 5.0) + ob_data = np.zeros((40, 40), dtype=np.float32) + ob_data[self._box] = 3.0 + led_data = self._extract(data, ob_data, np.zeros((40, 40), dtype=np.float32)) + assert led_data.ob_sum_color_value == pytest.approx(3.0 * ob_data[self._box].size) + + def test_negative_values_are_kept_in_the_sum(self): + # this is the whole point: the clipped sum cannot go below the pedestal it carries + data = np.zeros((40, 40)) + ob_data = np.full((40, 40), -2.0, dtype=np.float32) + led_data = self._extract(data, ob_data, np.zeros((40, 40), dtype=np.float32)) + assert led_data.ob_sum_color_value == pytest.approx(-2.0 * ob_data[self._box].size) + + def test_the_noise_adds_up_quadratically(self): + data = np.zeros((40, 40)) + noise = np.zeros((40, 40), dtype=np.float32) + noise[self._box] = 2.0 + led_data = self._extract(data, np.zeros((40, 40), dtype=np.float32), noise) + num_pixels = noise[self._box].size + assert led_data.ob_noise_sum_color_value == pytest.approx(2.0 * np.sqrt(num_pixels)) + + def test_pixels_outside_the_channel_carry_no_noise(self): + # the reader zeroes both arrays outside the channel, so a Bayer image sums only its own pixels + data = np.zeros((40, 40)) + noise = np.zeros((40, 40), dtype=np.float32) + noise[::2, ::2] = 2.0 + led_data = self._extract(data, np.zeros((40, 40), dtype=np.float32), noise) + num_channel_pixels = np.count_nonzero(noise[self._box]) + assert led_data.ob_noise_sum_color_value == pytest.approx(2.0 * np.sqrt(num_channel_pixels)) diff --git a/ledsa/tests/UnitTests/test_image_orientation.py b/ledsa/tests/UnitTests/test_image_orientation.py index 1a69f05..b276454 100644 --- a/ledsa/tests/UnitTests/test_image_orientation.py +++ b/ledsa/tests/UnitTests/test_image_orientation.py @@ -12,7 +12,8 @@ from ledsa.core import image_reading from ledsa.core.image_reading import (ORIENTATION_TRANSFORMATIONS, _apply_orientation, _orientation_from_rotation, - get_bayer_channel_mask, read_channel_data_and_mask_from_img, + get_bayer_channel_mask, + read_channel_data_mask_and_optical_black_from_img, read_channel_data_from_img, read_img_array_from_img) from ledsa.data_extraction.step_1_functions import find_search_areas @@ -118,7 +119,8 @@ def test_all_channels_are_read_upright(self, tmp_path): def test_channel_mask_fits_the_upright_array(self, tmp_path): stored = np.rot90(led_strip_img_array()) path = write_img(tmp_path / 'sideways.png', stored, orientation=SIDEWAYS) - channel_array, channel_mask = read_channel_data_and_mask_from_img(path, channel=0) + channel_array, channel_mask, _, _ = \ + read_channel_data_mask_and_optical_black_from_img(path, channel=0) assert channel_mask.shape == channel_array.shape and channel_mask.all() def test_rotation_option_overrides_the_metadata(self, tmp_path): @@ -193,7 +195,8 @@ def test_sideways_raw_is_read_upright(self, monkeypatch, stored): def test_bayer_pattern_is_turned_with_the_pixels(self, monkeypatch, stored): """The mask has to keep naming the color of every pixel, or the sums mix the channels.""" self.fake_imread(monkeypatch, stored, flip=6) - channel_array, channel_mask = read_channel_data_and_mask_from_img('img.CR3', channel=1) + channel_array, channel_mask, _, _ = \ + read_channel_data_mask_and_optical_black_from_img('img.CR3', channel=1) np.testing.assert_array_equal(channel_mask, np.rot90(get_bayer_channel_mask(self.FILTER_ARRAY, 1), -1)) np.testing.assert_array_equal(channel_array, np.where(channel_mask, np.rot90(stored, -1), 0)) diff --git a/ledsa/tests/UnitTests/test_sideband_background.py b/ledsa/tests/UnitTests/test_sideband_background.py index c0a4330..a11a698 100644 --- a/ledsa/tests/UnitTests/test_sideband_background.py +++ b/ledsa/tests/UnitTests/test_sideband_background.py @@ -286,9 +286,13 @@ def test_bayer_channel_mask(self, channel, colors): def test_raw_mask_keeps_channel_pixels_with_value_zero(self, monkeypatch): data = np.where(self.FILTER_ARRAY == 1, 0, 7) # the first green channel is clipped at 0 - monkeypatch.setattr(image_reading, '_read_img_array_from_raw_file', - lambda filename, channel, orientation=None: (data, self.FILTER_ARRAY)) - channel_array, mask = image_reading.read_channel_data_and_mask_from_img('img.CR3', 1) + # mocked at the one pass over the RAW file, so that every reader derived from it stays consistent. + # A black level of 0 and a white level of 255 leave the stored values as they are. + raw = image_reading.RawImage(data, self.FILTER_ARRAY, 0.0, 255, + image_reading.UPRIGHT_ORIENTATION, None, None) + monkeypatch.setattr(image_reading, '_read_raw_image', lambda filename, channel: raw) + channel_array, mask, _, _ = \ + image_reading.read_channel_data_mask_and_optical_black_from_img('img.CR3', 1) np.testing.assert_array_equal(mask, np.isin(self.FILTER_ARRAY, [1, 3])) np.testing.assert_array_equal(channel_array, np.where(mask, data, 0)) np.testing.assert_array_equal(channel_array, image_reading._read_channel_data_from_raw_file('img.CR3', 1)) @@ -297,11 +301,13 @@ def test_all_pixels_of_a_png_belong_to_the_channel(self, tmp_path): import matplotlib.pyplot as plt path = tmp_path / 'img.png' plt.imsave(path, np.random.default_rng(0).random((6, 8, 3))) - channel_array, mask = image_reading.read_channel_data_and_mask_from_img(str(path), 1) + channel_array, mask, ob_data, ob_noise = \ + image_reading.read_channel_data_mask_and_optical_black_from_img(str(path), 1) + assert ob_data is None and ob_noise is None # a PNG has no optical black np.testing.assert_array_equal(channel_array, image_reading.read_channel_data_from_img(str(path), 1)) assert mask.shape == channel_array.shape and mask.all() def test_unsupported_format_raises(self): from ledsa.core.image_integrity import CorruptedImageError with pytest.raises(CorruptedImageError, match='unsupported file extension'): - image_reading.read_channel_data_and_mask_from_img('img.xyz', 0) + image_reading.read_channel_data_mask_and_optical_black_from_img('img.xyz', 0) From 552b2fe7cf76df795cdfcb9e8e6b675712f311d6 Mon Sep 17 00:00:00 2001 From: Kristian Date: Wed, 28 Oct 2026 22:38:21 +0100 Subject: [PATCH 3/3] weight the LEDs of the linear solver by their signal to noise ratio With snr_full_weight the linear solver weights every LED with min(1, SNR / snr_full_weight). The standard deviation of the optical depth is about 1/SNR, so the weight is the inverse standard deviation, capped so that well measured LEDs keep the weight 1 and the meaning of lambda_reg does not change. LEDs at or below zero signal get the weight 0, and an image whose weights do not add up to one full LED yields NaN. Layers whose path length weighted mean LED weight stays below min_layer_weight become NaN as well: their coefficients would rest on LEDs near the noise level, where the linearised Beer-Lambert law is biased. The weighting needs the noise of the reference property, which only ob_sum_col_val has (NOISE_COLUMNS). snr_full_weight with any other reference property is refused with the columns that do have one, and the loader stops if the noise column is missing from the binary data. num_moving_avg_images smooths the reference property over N consecutive images with a centred moving average, and the noise becomes the one of that mean. It is the counterpart of average_images with n_summarize, which summarizes images into fewer rows and numbers them from 1 again: the moving average keeps every image with its own time and its own coefficient. Combining the two is refused, because the result would be smoothed and decimated at once. Unit tests pass, 337, with test_led_weighting.py for the weights, the discarded layers, the two refusals and the smoothing. Co-Authored-By: Claude Opus 5 --- ledsa/analysis/ConfigDataAnalysis.py | 23 ++- ledsa/analysis/ExperimentData.py | 19 +++ ledsa/analysis/ExtinctionCoefficients.py | 83 +++++++++- .../analysis/ExtinctionCoefficientsLinear.py | 46 +++++- ledsa/core/parser_arguments_run.py | 5 +- ledsa/tests/UnitTests/test_led_weighting.py | 151 ++++++++++++++++++ .../UnitTests/test_sideband_background.py | 3 + 7 files changed, 317 insertions(+), 13 deletions(-) create mode 100644 ledsa/tests/UnitTests/test_led_weighting.py diff --git a/ledsa/analysis/ConfigDataAnalysis.py b/ledsa/analysis/ConfigDataAnalysis.py index ef477fb..a53e126 100644 --- a/ledsa/analysis/ConfigDataAnalysis.py +++ b/ledsa/analysis/ConfigDataAnalysis.py @@ -58,7 +58,8 @@ def __init__(self, load_config_file=True, camera_position=None, num_layers=20, d reference_property='sum_col_val', average_images=False, n_summarize=2, solver='linear', weighting_preference=-6e-3, weighting_curvature=1e-2, - num_iterations=1000, lambda_reg=1e-3, reduction_files=None): + num_iterations=1000, lambda_reg=1e-3, snr_full_weight=None, min_layer_weight=0.5, + num_moving_avg_images=1, reduction_files=None): """ :param load_config_file: Determines whether to load the config file on initialization. Defaults to True. :type load_config_file: bool @@ -99,6 +100,18 @@ def __init__(self, load_config_file=True, camera_position=None, num_layers=20, d :type num_iterations: int :param lambda_reg: Regularization parameter for the linear solver. Defaults to 1e-3. :type lambda_reg: float + :param snr_full_weight: Linear solver only: signal-to-noise ratio from which on an LED gets the + full weight. Below it the weight falls linearly to zero at a signal of zero. Needs a reference + property with a noise column, which only ob_sum_col_val has. None weights all LEDs equally and + is the default. + :type snr_full_weight: float or None + :param min_layer_weight: Linear solver with snr_full_weight only: a layer whose path length + weighted mean LED weight stays below this value is set to NaN. Defaults to 0.5. + :type min_layer_weight: float + :param num_moving_avg_images: Linear solver only: number of consecutive images over which the + reference property is smoothed by a centred moving average, every image keeping its own time. + 1 leaves it as it is and is the default. Cannot be combined with average_images. + :type num_moving_avg_images: int :param reduction_files: Path of the reduction factor file of the temperature correction per camera channel. Channels without an entry are analysed without the correction, which is the default for every channel. @@ -133,6 +146,14 @@ def __init__(self, load_config_file=True, camera_position=None, num_layers=20, d self['DEFAULT'][' num_iterations'] = str(num_iterations) self.set('DEFAULT', ' # Options for linear solver (ignored when nonlinear solver is used) ') self['DEFAULT'][' lambda_reg'] = str(lambda_reg) + self.set('DEFAULT', ' # SNR from which on an LED gets the full weight, None for equal weights.') + self.set('DEFAULT', ' # Needs a reference property with a noise column (ob_sum_col_val).') + self['DEFAULT'][' snr_full_weight'] = str(snr_full_weight) + self.set('DEFAULT', ' # Layers whose mean LED weight stays below this value become NaN') + self['DEFAULT'][' min_layer_weight'] = str(min_layer_weight) + self.set('DEFAULT', ' # Images for a centred moving average of the reference property, 1 is off.') + self.set('DEFAULT', ' # Every image keeps its own time, unlike average_images with n_summarize.') + self['DEFAULT'][' num_moving_avg_images'] = str(num_moving_avg_images) self['experiment_geometry'] = {} self.set('experiment_geometry', '# Global X Y Z position of the camera ') diff --git a/ledsa/analysis/ExperimentData.py b/ledsa/analysis/ExperimentData.py index 05675ab..2fa32d9 100644 --- a/ledsa/analysis/ExperimentData.py +++ b/ledsa/analysis/ExperimentData.py @@ -36,6 +36,15 @@ class ExperimentData: :type num_ref_images: int :ivar ref_img_indices: Indices of reference images to use. If None, use num_ref_imgs. :type ref_img_indices: list[int] or None + :ivar snr_full_weight: Signal-to-noise ratio from which on an LED gets the full weight in the linear + solver. If None, all LEDs are weighted equally. + :type snr_full_weight: float or None + :ivar min_layer_weight: Only with snr_full_weight: a layer whose path length weighted mean LED weight + stays below this value is set to NaN. + :type min_layer_weight: float + :ivar num_moving_avg_images: Number of consecutive images for a centred moving average of the reference + property. 1 leaves it as it is. + :type num_moving_avg_images: int :ivar average_images: Intensities are computed as an average from consecutive images. :type average_images: bool :ivar n_summarize: Number of consecutive images summarized into one when averaging. @@ -66,6 +75,9 @@ def __init__(self, load_config_file=True): self.average_images = False self.n_summarize = 2 self.lambda_reg = None + self.snr_full_weight = None + self.min_layer_weight = 0.5 + self.num_moving_avg_images = 1 self.reference_properties = None self.reduction_files = {} self.merge_led_arrays = None @@ -94,6 +106,13 @@ def load_config_parameters(self) -> None: self.num_iterations = int(config_analysis['DEFAULT']['num_iterations']) elif self.solver == 'linear': self.lambda_reg = float(config_analysis['DEFAULT']['lambda_reg']) + if config_analysis.has_option('DEFAULT', 'snr_full_weight') and \ + config_analysis['DEFAULT']['snr_full_weight'].strip().lower() != 'none': + self.snr_full_weight = float(config_analysis['DEFAULT']['snr_full_weight']) + if config_analysis.has_option('DEFAULT', 'min_layer_weight'): + self.min_layer_weight = float(config_analysis['DEFAULT']['min_layer_weight']) + if config_analysis.has_option('DEFAULT', 'num_moving_avg_images'): + self.num_moving_avg_images = int(config_analysis['DEFAULT']['num_moving_avg_images']) self.reference_properties = config_analysis.get_reference_properties() self.reduction_files = {channel: config_analysis.get_reduction_file(channel) for channel in (self.channels or [])} diff --git a/ledsa/analysis/ExtinctionCoefficients.py b/ledsa/analysis/ExtinctionCoefficients.py index b982424..b402de9 100644 --- a/ledsa/analysis/ExtinctionCoefficients.py +++ b/ledsa/analysis/ExtinctionCoefficients.py @@ -8,7 +8,7 @@ from ledsa.analysis.Experiment import Experiment, Layers, Camera from ledsa.analysis.temperature_correction import TemperatureCorrection -from ledsa.core.file_handling import (NUM_IMGS_PER_AVERAGE, create_analysis_infos_avg, extend_hdf, +from ledsa.core.file_handling import (NOISE_COLUMNS, NUM_IMGS_PER_AVERAGE, create_analysis_infos_avg, extend_hdf, get_extinction_coefficients_file_path, read_experiment_times, read_hdf, read_hdf_avg) from ledsa.data_extraction.data_integrity import check_intensity_normalization @@ -53,7 +53,7 @@ class ExtinctionCoefficients(ABC): def __init__(self, experiment, reference_property='sum_col_val', num_ref_imgs=10, ref_img_indices=None, average_images=False, n_summarize=NUM_IMGS_PER_AVERAGE, - temperature_correction=None): + temperature_correction=None, snr_full_weight=None, num_moving_avg_images=1): """ :param experiment: Object representing the experimental setup. :type experiment: Experiment @@ -70,7 +70,23 @@ def __init__(self, experiment, reference_property='sum_col_val', num_ref_imgs=10 :param temperature_correction: Reduction factors of the LED intensities of the channel. None leaves the intensities as they were extracted. :type temperature_correction: TemperatureCorrection or None + :param snr_full_weight: Signal-to-noise ratio from which on an LED gets the full weight in the + solver. Below it the weight falls linearly to zero at a signal of zero. Needs a reference + property with a noise column. If None, all LEDs are weighted equally. + :type snr_full_weight: float or None + :param num_moving_avg_images: Number of consecutive images over which the reference property is + smoothed by a centred moving average before the normalisation. 1 leaves it as it is. Unlike + average_images this keeps every image with its own time and does not renumber the images. + :type num_moving_avg_images: int + :raises ValueError: If snr_full_weight is given for a reference property without a noise column, or + if average_images and num_moving_avg_images are combined. """ + if average_images and int(num_moving_avg_images) > 1: + raise ValueError('average_images and num_moving_avg_images cannot be combined. ' + 'average_images summarizes n_summarize images into one and numbers the ' + 'result from 1 again, while num_moving_avg_images smooths the reference ' + 'property over neighbouring images and keeps every image with its own ' + 'time. Pick the one that suits the time resolution you need.') self.coefficients_per_image_and_layer = [] self.experiment = experiment self.reference_property = reference_property @@ -87,10 +103,19 @@ def __init__(self, experiment, reference_property='sum_col_val', num_ref_imgs=10 self.temperature_corrected = False self.solver = None self.lambda_reg = None + self.snr_full_weight = snr_full_weight + self.num_moving_avg_images = int(num_moving_avg_images) + self.noise_property = None + if snr_full_weight is not None: + if reference_property not in NOISE_COLUMNS: + raise ValueError(f'snr_full_weight needs a reference property with a noise column, but ' + f'{reference_property} has none. Columns with one: ' + f'{", ".join(sorted(NOISE_COLUMNS))}.') + self.noise_property = NOISE_COLUMNS[reference_property] def __str__(self): out = str(self.experiment) + \ - f'Reference_property: {self.reference_property}, Ref_img_indices: {self.ref_img_indices}, Solver: {self.solver}, Lambda_reg: {self.lambda_reg}, Temperature_correction: {self.temperature_correction_info()}, LEDSA {version("ledsa")}\n' + f'Reference_property: {self.reference_property}, Ref_img_indices: {self.ref_img_indices}, Snr_full_weight: {self.snr_full_weight}, Min_layer_weight: {getattr(self, "min_layer_weight", None) if self.snr_full_weight is not None else None}, Num_moving_avg_images: {self.num_moving_avg_images}, Solver: {self.solver}, Lambda_reg: {self.lambda_reg}, Temperature_correction: {self.temperature_correction_info()}, LEDSA {version("ledsa")}\n' return out def temperature_correction_info(self) -> str: @@ -117,7 +142,11 @@ def calc_and_set_coefficients(self) -> None: rel_intensities = single_img_array / self.ref_intensities # Calculate the extinction coefficients depending on child class used - sigmas = self.calc_coefficients_of_img(rel_intensities) + if self.noise_property is None: + sigmas = self.calc_coefficients_of_img(rel_intensities) + else: + snr = single_img_array / single_img_data[self.noise_property].to_numpy() + sigmas = self.calc_coefficients_of_img(rel_intensities, snr) self.coefficients_per_image_and_layer.append(sigmas) def calc_and_set_coefficients_mp(self, cores=4) -> None: @@ -134,7 +163,12 @@ def calc_and_set_coefficients_mp(self, cores=4) -> None: # Calculate the extinction coefficients depending on child class used pool = Pool(processes=cores) - sigmas = pool.map(self.calc_coefficients_of_img, rel_intensities) + if self.noise_property is None: + sigmas = pool.map(self.calc_coefficients_of_img, rel_intensities) + else: + snr = img_property_array / multiindex_series_to_nparray( + self.calculated_img_data[self.noise_property]) + sigmas = pool.starmap(self.calc_coefficients_of_img, zip(rel_intensities, snr)) pool.close() self.coefficients_per_image_and_layer = sigmas @@ -192,7 +226,14 @@ def load_img_data(self) -> None: raise ValueError(f'{self.reference_property} is not among the columns of the binary data of ' f'channel {self.experiment.channel} ' f'({", ".join(map(str, img_data.columns))}). {hint}') - img_data_cropped = img_data[['led_array_id', self.reference_property]] + columns = ['led_array_id', self.reference_property] + if self.noise_property is not None: + if self.noise_property not in img_data.columns: + raise ValueError(f'{self.noise_property} is needed for snr_full_weight but is not among ' + f'the columns of the binary data of channel {self.experiment.channel}. ' + f'Run step 3 again with a LEDSA version that writes it.') + columns.append(self.noise_property) + img_data_cropped = img_data[columns] self.calculated_img_data = img_data_cropped[img_data_cropped['led_array_id'] == self.experiment.led_array] if self.calculated_img_data.empty: exit(f"Apparently there are no intensity values for led array {self.experiment.led_array}!") @@ -210,6 +251,36 @@ def load_img_data(self) -> None: # gapless range. They are kept to pair every computed coefficient with the time of its own image. self.analysed_img_ids = sorted( self.calculated_img_data.index.get_level_values('img_id').unique().tolist()) + if self.num_moving_avg_images > 1: + self.apply_moving_average() + + def apply_moving_average(self) -> None: + """ + Replace the reference property of every LED by its centred moving average over + num_moving_avg_images consecutive analysed images, fewer at the start and the end of the series. + The noise, where it is used, becomes the standard deviation of that mean, so the root of the summed + variances divided by the number of images. + + Every image keeps its own time and its own row, so the series is smoothed rather than decimated: + the time resolution becomes num_moving_avg_images times the image interval while the number of + extinction coefficients stays the same. average_images with n_summarize is the other way round. + """ + n = self.num_moving_avg_images + columns = [self.reference_property] + ([self.noise_property] if self.noise_property else []) + averaged = {} + for column in columns: + wide = self.calculated_img_data[column].unstack('led_id') + if column == self.noise_property: + count = wide.notna().astype(float).rolling(n, center=True, min_periods=1).sum() + wide = np.sqrt((wide ** 2).rolling(n, center=True, min_periods=1).sum()) / count + else: + wide = wide.rolling(n, center=True, min_periods=1).mean() + averaged[column] = wide.stack() + data = self.calculated_img_data.copy() + for column, values in averaged.items(): + data[column] = values.reindex(data.index).to_numpy() + self.calculated_img_data = data + print(f'Reference property smoothed over {n} consecutive images (centred moving average).') def apply_temperature_correction(self) -> None: """ diff --git a/ledsa/analysis/ExtinctionCoefficientsLinear.py b/ledsa/analysis/ExtinctionCoefficientsLinear.py index dfd34ea..fb63535 100644 --- a/ledsa/analysis/ExtinctionCoefficientsLinear.py +++ b/ledsa/analysis/ExtinctionCoefficientsLinear.py @@ -19,7 +19,8 @@ class ExtinctionCoefficientsLinear(ExtinctionCoefficients): :vartype solver: str """ - def __init__(self, experiment, reference_property='sum_col_val', num_ref_imgs=10, average_images=False, ref_img_indices=None, lambda_reg=1e-3, n_summarize=2, temperature_correction=None): + def __init__(self, experiment, reference_property='sum_col_val', num_ref_imgs=10, average_images=False, ref_img_indices=None, lambda_reg=1e-3, n_summarize=2, temperature_correction=None, + snr_full_weight=None, min_layer_weight=0.5, num_moving_avg_images=1): """ Initialize the ExtinctionCoefficientsLinear object. @@ -40,20 +41,41 @@ def __init__(self, experiment, reference_property='sum_col_val', num_ref_imgs=10 :param temperature_correction: Reduction factors of the LED intensities of the channel. None leaves the intensities as they were extracted. :type temperature_correction: TemperatureCorrection or None + :param snr_full_weight: Signal-to-noise ratio from which on an LED gets the full weight. Below it + the weight falls linearly to zero at a signal of zero. If None, all LEDs are weighted equally. + :type snr_full_weight: float or None + :param min_layer_weight: Only with snr_full_weight: a layer is set to NaN if the mean weight of the + LEDs crossing it, weighted with their path length in the layer, stays below this value. + :type min_layer_weight: float + :param num_moving_avg_images: Number of consecutive images for a centred moving average of the + reference property. + :type num_moving_avg_images: int """ + self.min_layer_weight = min_layer_weight super().__init__(experiment, reference_property, num_ref_imgs, ref_img_indices, average_images, - n_summarize, temperature_correction) + n_summarize, temperature_correction, snr_full_weight, num_moving_avg_images) self.lambda_reg = lambda_reg self.solver = 'linear' - def calc_coefficients_of_img(self, rel_intensities: np.ndarray) -> np.ndarray: + def calc_coefficients_of_img(self, rel_intensities: np.ndarray, snr: np.ndarray = None) -> np.ndarray: """ Calculate the extinction coefficients for a single image using a linearized approach. This method solves the linear system derived from the Beer-Lambert law using non-negative least squares with Tikhonov regularization. + If snr is given, every LED is weighted with min(1, snr / snr_full_weight). The standard deviation + of the optical depth -ln(I_e/I_0) is about 1/snr, so the weight is the inverse standard deviation, + capped so that well measured LEDs keep the weight 1 and the meaning of lambda_reg does not change. + LEDs with a signal at or below zero get the weight 0. If the weights add up to less than one LED + with full weight, no coefficients can be computed and NaN is returned. Layers in which the path + length weighted mean weight of the crossing LEDs stays below min_layer_weight are set to NaN as + well: their coefficients would rest on LEDs near the noise level, where the linearised + Beer-Lambert law is biased. + :param rel_intensities: Array of relative (normalized) LED intensities (I_e/I_0). :type rel_intensities: np.ndarray + :param snr: Signal-to-noise ratio of every LED, or None for equal weights. + :type snr: np.ndarray or None :return: Array of the computed extinction coefficients (sigmas). :rtype: np.ndarray """ @@ -68,6 +90,15 @@ def calc_coefficients_of_img(self, rel_intensities: np.ndarray) -> np.ndarray: # Get dimensions of the problem n_leds, n_layers = self.distances_per_led_and_layer.shape + # Weights of the LEDs: the inverse standard deviation of the optical depth, capped at 1 + if snr is None: + weights = np.ones(n_leds) + else: + weights = np.clip(np.nan_to_num(np.asarray(snr, dtype=float) / self.snr_full_weight, nan=0.0), + 0.0, 1.0) + if weights.sum() < 1.0: + return np.full(n_layers, np.nan) + # Build finite-difference matrix L for second derivative (curvature penalty) # This matrix approximates the second derivative and is used for regularization # to enforce smoothness in the solution @@ -86,15 +117,20 @@ def calc_coefficients_of_img(self, rel_intensities: np.ndarray) -> np.ndarray: # Augment the system with regularization # A_aug = [A; sqrt(lambda)*L], where A is the distance matrix - A_aug = np.vstack((self.distances_per_led_and_layer, sqrt_lambda * L)) + A_aug = np.vstack((self.distances_per_led_and_layer * weights[:, None], sqrt_lambda * L)) # Augment the right-hand side: b_aug = [b; 0] - b_aug = np.hstack((optical_depths, np.zeros(L.shape[0]))) + b_aug = np.hstack((optical_depths * weights, np.zeros(L.shape[0]))) # Solve the non-negative least squares problem: min ||A_aug*x - b_aug||^2 subject to x >= 0 sigmas, residuals = nnls(A_aug, b_aug) # sigmas, residuals, rank, s = np.linalg.lstsq(A_aug, b_aug, rcond=None) + if snr is not None and self.min_layer_weight is not None: + path = self.distances_per_led_and_layer + layer_weight = (weights[:, None] * path).sum(axis=0) / np.maximum(path.sum(axis=0), 1e-12) + sigmas = np.where(layer_weight < self.min_layer_weight, np.nan, sigmas) + return sigmas # For investigation diff --git a/ledsa/core/parser_arguments_run.py b/ledsa/core/parser_arguments_run.py index 84f30b3..8e5fa95 100644 --- a/ledsa/core/parser_arguments_run.py +++ b/ledsa/core/parser_arguments_run.py @@ -285,7 +285,10 @@ def extionction_coefficient_calculation(args) -> None: average_images=ex_data.average_images, n_summarize=ex_data.n_summarize, temperature_correction=temperature_corrections.get(channel), - lambda_reg=ex_data.lambda_reg) + lambda_reg=ex_data.lambda_reg, + snr_full_weight=ex_data.snr_full_weight, + min_layer_weight=ex_data.min_layer_weight, + num_moving_avg_images=ex_data.num_moving_avg_images) else: raise ValueError(f"Invalid solver type '{solver}'. Must be 'linear' or 'nonlinear'.") if ex_data.n_cpus > 1: diff --git a/ledsa/tests/UnitTests/test_led_weighting.py b/ledsa/tests/UnitTests/test_led_weighting.py new file mode 100644 index 0000000..a64df74 --- /dev/null +++ b/ledsa/tests/UnitTests/test_led_weighting.py @@ -0,0 +1,151 @@ +from types import SimpleNamespace + +import numpy as np +import pandas as pd +import pytest + +from ledsa.analysis.ExtinctionCoefficientsLinear import ExtinctionCoefficientsLinear + +# 4 LEDs over 3 layers: the first three cross one layer each, the fourth crosses all of them +DISTANCES = np.array([[1.0, 0.0, 0.0], + [0.0, 1.0, 0.0], + [0.0, 0.0, 1.0], + [1.0, 1.0, 1.0]]) +SIGMAS = np.array([0.1, 0.2, 0.3]) +REL_INTENSITIES = np.exp(-DISTANCES @ SIGMAS) +SNR_FULL_WEIGHT = 10.0 + + +def make_solver(min_layer_weight=None) -> ExtinctionCoefficientsLinear: + """A linear solver with the geometry above, built without __init__ as the config is not needed.""" + solver = ExtinctionCoefficientsLinear.__new__(ExtinctionCoefficientsLinear) + solver.distances_per_led_and_layer = DISTANCES + solver.lambda_reg = 1e-6 + solver.snr_full_weight = SNR_FULL_WEIGHT + solver.min_layer_weight = min_layer_weight + return solver + + +class TestWhichReferencePropertiesCanBeWeighted: + def test_a_property_without_a_noise_column_is_refused(self): + experiment = SimpleNamespace(channel=0, path='.', led_array=0) + with pytest.raises(ValueError, match=r'snr_full_weight needs a reference property with a noise ' + r'column.*sidebg_sum_col_val has none.*ob_sum_col_val'): + ExtinctionCoefficientsLinear(experiment, reference_property='sidebg_sum_col_val', + snr_full_weight=SNR_FULL_WEIGHT) + + def test_the_unclipped_sum_brings_its_noise_column(self): + experiment = SimpleNamespace(channel=0, path='.', led_array=0) + coefficients = ExtinctionCoefficientsLinear(experiment, reference_property='ob_sum_col_val', + snr_full_weight=SNR_FULL_WEIGHT) + assert coefficients.noise_property == 'ob_noise_sum_col_val' + + def test_without_the_option_no_noise_column_is_needed(self): + experiment = SimpleNamespace(channel=0, path='.', led_array=0) + assert ExtinctionCoefficientsLinear(experiment, reference_property='sum_col_val').noise_property is None + + +class TestSmoothingVersusAveraging: + def test_the_two_ways_of_averaging_cannot_be_combined(self): + # one summarizes images into fewer rows and renumbers them, the other smooths and keeps the rows + experiment = SimpleNamespace(channel=0, path='.', led_array=0) + with pytest.raises(ValueError, match=r'average_images and num_moving_avg_images cannot be ' + r'combined'): + ExtinctionCoefficientsLinear(experiment, average_images=True, num_moving_avg_images=3) + + def test_smoothing_alone_is_allowed(self): + experiment = SimpleNamespace(channel=0, path='.', led_array=0) + assert ExtinctionCoefficientsLinear(experiment, num_moving_avg_images=3).num_moving_avg_images == 3 + + +class TestLedWeights: + def test_equal_weights_reproduce_the_unweighted_solution(self): + solver = make_solver() + snr = np.full(DISTANCES.shape[0], 5 * SNR_FULL_WEIGHT) # every LED well above the full weight + np.testing.assert_allclose(solver.calc_coefficients_of_img(REL_INTENSITIES, snr), + solver.calc_coefficients_of_img(REL_INTENSITIES), atol=1e-9) + + def test_the_unweighted_solution_recovers_the_coefficients(self): + np.testing.assert_allclose(make_solver().calc_coefficients_of_img(REL_INTENSITIES), SIGMAS, + atol=1e-3) + + def test_an_led_without_signal_does_not_influence_the_result(self): + # its weight is 0, so whatever it reports must not reach the coefficients + solver = make_solver() + snr = np.array([0.0, 5 * SNR_FULL_WEIGHT, 5 * SNR_FULL_WEIGHT, 5 * SNR_FULL_WEIGHT]) + first, second = REL_INTENSITIES.copy(), REL_INTENSITIES.copy() + first[0], second[0] = 0.9, 0.01 + np.testing.assert_allclose(solver.calc_coefficients_of_img(first, snr), + solver.calc_coefficients_of_img(second, snr), atol=1e-9) + + def test_a_nan_signal_to_noise_ratio_becomes_weight_zero(self): + solver = make_solver() + snr = np.array([np.nan, 5 * SNR_FULL_WEIGHT, 5 * SNR_FULL_WEIGHT, 5 * SNR_FULL_WEIGHT]) + assert np.isfinite(solver.calc_coefficients_of_img(REL_INTENSITIES, snr)).all() + + def test_too_little_signal_in_the_whole_image_gives_nan(self): + # the weights add up to less than one LED with full weight, so nothing can be computed + solver = make_solver() + snr = np.full(DISTANCES.shape[0], 0.1 * SNR_FULL_WEIGHT) + assert np.isnan(solver.calc_coefficients_of_img(REL_INTENSITIES, snr)).all() + + +class TestLayersRestingOnNoise: + SNR = np.array([0.0, 0.0, 0.0, 2 * SNR_FULL_WEIGHT]) # only the LED crossing every layer is measured + + def test_layers_below_the_limit_become_nan(self): + # every layer carries the path weighted mean weight 0.5, from the one LED that crosses all of them + assert np.isnan(make_solver(min_layer_weight=0.6).calc_coefficients_of_img(REL_INTENSITIES, + self.SNR)).all() + + def test_layers_above_the_limit_are_kept(self): + assert np.isfinite(make_solver(min_layer_weight=0.4).calc_coefficients_of_img(REL_INTENSITIES, + self.SNR)).all() + + def test_without_the_limit_nothing_is_discarded(self): + assert np.isfinite(make_solver().calc_coefficients_of_img(REL_INTENSITIES, self.SNR)).all() + + +class TestMovingAverage: + NUM_IMAGES, NUM_LEDS = 5, 2 + + def make_coefficients(self, values, noise=None): + coefficients = ExtinctionCoefficientsLinear.__new__(ExtinctionCoefficientsLinear) + coefficients.reference_property = 'ob_sum_col_val' + coefficients.noise_property = 'ob_noise_sum_col_val' if noise is not None else None + coefficients.num_moving_avg_images = 3 + index = pd.MultiIndex.from_product([range(1, self.NUM_IMAGES + 1), range(self.NUM_LEDS)], + names=['img_id', 'led_id']) + data = {'led_array_id': 0, 'ob_sum_col_val': values} + if noise is not None: + data['ob_noise_sum_col_val'] = noise + coefficients.calculated_img_data = pd.DataFrame(data, index=index) + return coefficients + + def test_every_image_keeps_its_own_row(self): + coefficients = self.make_coefficients(np.arange(self.NUM_IMAGES * self.NUM_LEDS, dtype=float)) + coefficients.apply_moving_average() + assert len(coefficients.calculated_img_data) == self.NUM_IMAGES * self.NUM_LEDS + assert coefficients.calculated_img_data.index.get_level_values('img_id').nunique() == self.NUM_IMAGES + + def test_a_spike_is_smoothed_over_its_neighbours(self): + values = np.zeros(self.NUM_IMAGES * self.NUM_LEDS) + values[2 * self.NUM_LEDS] = 30.0 # led 0 of the middle image + coefficients = self.make_coefficients(values) + coefficients.apply_moving_average() + smoothed = coefficients.calculated_img_data['ob_sum_col_val'].unstack('led_id')[0].to_numpy() + np.testing.assert_allclose(smoothed, [0.0, 10.0, 10.0, 10.0, 0.0]) + + def test_the_ends_average_over_fewer_images(self): + coefficients = self.make_coefficients(np.ones(self.NUM_IMAGES * self.NUM_LEDS)) + coefficients.apply_moving_average() + assert coefficients.calculated_img_data['ob_sum_col_val'].to_numpy() == pytest.approx(1.0) + + def test_the_noise_becomes_the_one_of_the_mean(self): + coefficients = self.make_coefficients(np.ones(self.NUM_IMAGES * self.NUM_LEDS), + noise=np.full(self.NUM_IMAGES * self.NUM_LEDS, 3.0)) + coefficients.apply_moving_average() + noise = coefficients.calculated_img_data['ob_noise_sum_col_val'].unstack('led_id')[0].to_numpy() + # three images in the middle, two at the ends: sqrt(n * sigma^2) / n = sigma / sqrt(n) + np.testing.assert_allclose(noise, [3 / np.sqrt(2), 3 / np.sqrt(3), 3 / np.sqrt(3), + 3 / np.sqrt(3), 3 / np.sqrt(2)]) diff --git a/ledsa/tests/UnitTests/test_sideband_background.py b/ledsa/tests/UnitTests/test_sideband_background.py index a11a698..1c9ccd5 100644 --- a/ledsa/tests/UnitTests/test_sideband_background.py +++ b/ledsa/tests/UnitTests/test_sideband_background.py @@ -261,6 +261,9 @@ def coefficients(self, monkeypatch): coefficients = ExtinctionCoefficientsLinear.__new__(ExtinctionCoefficientsLinear) coefficients.experiment = SimpleNamespace(channel=0, path='.', led_array=0) coefficients.average_images = False + # set in __init__: the property analysed here has no noise column and is not smoothed + coefficients.noise_property = None + coefficients.num_moving_avg_images = 1 return coefficients def test_nan_values_stop_the_analysis(self, coefficients):