diff --git a/README.md b/README.md index 61d0b66..25a899c 100644 --- a/README.md +++ b/README.md @@ -120,6 +120,37 @@ defines the LED search areas, the first image of the experiment (`first_img_experiment_id`), which defines the experiment start time, and the time reference image. +### Background correction (optional) + +By default the extinction coefficients are computed from `sum_col_val`, the plain sum of +the pixel values in the search area of each LED. Step 3 additionally writes +background-corrected sums, which can be selected with `reference_property` in +`config_analysis.ini`: + +| `reference_property` | Background estimate | +|---|---| +| `sum_col_val` (default) | none | +| `bgsub_sum_col_val` | median of the pixels > 0 in a 2 px border ring of the search area, subtracted pixel-wise | +| `bgsub2_sum_col_val` | mean over the channel pixels (including 0) of the same border ring, subtracted as N·b | +| `sidebg_sum_col_val` | mean over the channel pixels in side bands next to the LED strip, subtracted as N·B | + +N is the number of channel pixels in the search area. For small search areas the border +ring lies in the halo of the LED itself, whereas the side bands sample the additive +pedestal next to the LED: stray light, light from the smoke and noise clipped at 0. Under +dense smoke the box sum consists mostly of this pedestal. + +The side bands run parallel to the strip, 25 to 60 px across it (`sideband_distances` in +`config.ini`); pixels close to other search areas are excluded. The direction of the strip +is estimated for each LED from its nearest neighbours in the same LED array. +`sidebg_pedestal` holds the subtracted N·B. `sidebg_sum_col_val` is NaN if the direction +cannot be estimated (an LED array with a single LED) or if the side bands contain no pixel, +for example because the image ends closer to the strip. The analysis stops with an error +naming these LEDs if the selected `reference_property` contains NaN. + +`bgsub2_sum_col_val` and `sidebg_sum_col_val` are kept at least at `sideband_floor` +(default 5 %) of `sum_col_val`; values at this limit carry no information on the +extinction. + ## Demo LEDSA includes a demo that demonstrates its functionality using sample data. diff --git a/doc/cli_documentation.rst b/doc/cli_documentation.rst index fe73ae6..0533065 100644 --- a/doc/cli_documentation.rst +++ b/doc/cli_documentation.rst @@ -34,7 +34,7 @@ Data Extraction - Step 3: Analyse all images to extract pixel values inside the search areas and fit them to an LED model (SLOW!). - ``--r`` (Use the red channel, default), ``--g`` (Use the green channel), ``--b`` (Use the blue channel), ``-rgb`` (Run for each channel) * - ``-s3_fast``, ``--step_3_fast`, ``--analyse_photo_fast`` - - Step 3: Analyse all images to extract pixel values and calculate the accumulated values inside the search areas (FAST!). + - Step 3: Analyse all images to extract pixel values and calculate the accumulated values inside the search areas (FAST!). Besides the plain sum, the background-corrected sums ``bgsub_sum_col_val``, ``bgsub2_sum_col_val`` and ``sidebg_sum_col_val`` are written (see ``reference_property`` in the analysis configuration). - -- * - ``-re``, ``--restart`` - Restart step 3 if it was previously interrupted. Only the images that have not been analysed are taken into account. diff --git a/doc/data_extraction.rst b/doc/data_extraction.rst index b1416d1..2dae0da 100644 --- a/doc/data_extraction.rst +++ b/doc/data_extraction.rst @@ -7,6 +7,14 @@ Data Extraction :show-inheritance: .. automodule:: ledsa.data_extraction.LEDAnalysisData + :members: + :undoc-members: + :show-inheritance: + +Side Band Pedestal +------------------ + +.. automodule:: ledsa.data_extraction.sideband_background :members: :undoc-members: :show-inheritance: \ No newline at end of file diff --git a/ledsa/analysis/ExtinctionCoefficients.py b/ledsa/analysis/ExtinctionCoefficients.py index 4301c62..f1a9dd0 100644 --- a/ledsa/analysis/ExtinctionCoefficients.py +++ b/ledsa/analysis/ExtinctionCoefficients.py @@ -128,6 +128,7 @@ def load_img_data(self) -> None: """ Load processed image data from binary file + :raises ValueError: If the reference property contains NaN values. """ if self.average_images: img_data = read_hdf_avg(self.experiment.channel, path=self.experiment.path) @@ -138,6 +139,16 @@ def load_img_data(self) -> None: 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}!") + # The solvers cannot handle NaN: the linear one stops with an unspecific error, the nonlinear one returns + # zeros without a warning. sidebg_sum_col_val is NaN for LEDs without side band pixels, for example. + nan_values = self.calculated_img_data[self.reference_property].isna() + if nan_values.any(): + nan_led_ids = sorted(self.calculated_img_data[nan_values].index.get_level_values('led_id').unique()) + led_ids_text = ', '.join(map(str, nan_led_ids[:10])) + (', ...' if len(nan_led_ids) > 10 else '') + raise ValueError(f'{int(nan_values.sum())} values of {self.reference_property} of led array ' + f'{self.experiment.led_array} are NaN (led ids {led_ids_text}). For sidebg_sum_col_val ' + f'this means that no side band pixels were found for these LEDs, see sideband_distances ' + f'in config.ini.') # Images without analysis results are missing from the binary data, so the IDs are not necessarily a # gapless range. They are kept to pair every computed coefficient with the time of its own image. self.analysed_img_ids = sorted( diff --git a/ledsa/core/ConfigData.py b/ledsa/core/ConfigData.py index 951a4f0..741c8d7 100644 --- a/ledsa/core/ConfigData.py +++ b/ledsa/core/ConfigData.py @@ -17,7 +17,8 @@ def __init__(self, load_config_file=True, img_directory=None, search_area_radius time_diff_to_image_time=None, img_name_string=None, num_img_overflow=None, first_img_experiment_id=None, last_img_experiment_id=None, ref_img_id=None, ignore_led_indices=None, led_array_edge_indices=None, led_array_edge_coordinates=None, first_img_analysis_id=None, - last_img_analysis_id=None, num_skip_imgs=0, num_skip_leds=0, merge_led_array_indices=None): # TODO: merge LED arrays + last_img_analysis_id=None, num_skip_imgs=0, num_skip_leds=0, merge_led_array_indices=None, + sideband_distances='25 60', sideband_floor=0.05, sideband_exclusion_margin=None): # TODO: merge LED arrays """ :param load_config_file: Determines whether to load the config file on initialization. Defaults to True. :type load_config_file: bool @@ -76,6 +77,15 @@ def __init__(self, load_config_file=True, img_directory=None, search_area_radius :type num_skip_leds: int :param merge_led_array_indices: Flag to merge LED arrays for analysis. Defaults to None. TODO: not a flag but list of arrays to merge? :type merge_led_array_indices: bool or None + :param sideband_distances: Inner and outer distance in pixels across the LED strip of the side bands used to + estimate the pedestal for sidebg_sum_col_val. Defaults to '25 60'. + :type sideband_distances: str + :param sideband_floor: Lower limit of sidebg_sum_col_val and bgsub2_sum_col_val as fraction of sum_col_val. + Defaults to 0.05. + :type sideband_floor: float + :param sideband_exclusion_margin: Side band pixels within search_area_radius + margin of any other search area + are excluded. Defaults to None, i.e. margin = search_area_radius. + :type sideband_exclusion_margin: int or None """ cp.ConfigParser.__init__(self, allow_no_value=True) if load_config_file: @@ -150,6 +160,13 @@ def __init__(self, load_config_file=True, img_directory=None, search_area_radius self['analyse_photo'][' num_skip_imgs'] = str(num_skip_imgs) self.set('analyse_photo', ' # In the analysis, only LEDs are considered with id dividable by skip_leds + 1.') self['analyse_photo'][' num_skip_leds'] = str(num_skip_leds) + self.set('analyse_photo', ' # Pedestal estimation in side bands parallel to the LED strip (sidebg_sum_col_val):') + self.set('analyse_photo', ' # inner and outer distance of the side bands across the strip in pixels,') + self['analyse_photo'][' sideband_distances'] = str(sideband_distances) + self.set('analyse_photo', ' # lower limit of the corrected sum as fraction of the plain sum,') + self['analyse_photo'][' sideband_floor'] = str(sideband_floor) + self.set('analyse_photo', ' # pixels within search_area_radius + margin of other search areas are excluded, None: margin = radius.') + self['analyse_photo'][' sideband_exclusion_margin'] = str(sideband_exclusion_margin) with open('config.ini', 'w', encoding='utf-8') as configfile: self.write(configfile) diff --git a/ledsa/core/file_handling.py b/ledsa/core/file_handling.py index 0defc07..25812bb 100644 --- a/ledsa/core/file_handling.py +++ b/ledsa/core/file_handling.py @@ -289,15 +289,19 @@ def create_binary_data(channel: int, path: str = '.') -> None: # prepended when they are loaded, so a file with n columns has n + 1 names. 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"] 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"] # known layouts by the number of names (including img_id) COLUMN_LAYOUTS = { len(BASE_COLUMNS): BASE_COLUMNS, # -s3_fast, before the local background subtraction - len(BASE_COLUMNS + BGSUB_COLUMNS): BASE_COLUMNS + BGSUB_COLUMNS, # -s3_fast + len(BASE_COLUMNS + BGSUB_COLUMNS): BASE_COLUMNS + BGSUB_COLUMNS, # -s3_fast, before the side band pedestal + len(BASE_COLUMNS + BGSUB_COLUMNS + SIDEBG_COLUMNS): BASE_COLUMNS + BGSUB_COLUMNS + SIDEBG_COLUMNS, # -s3_fast len(BASE_COLUMNS + BGSUB_COLUMNS + LED_CENTER_COLUMNS + FIT_COLUMNS): - BASE_COLUMNS + BGSUB_COLUMNS + LED_CENTER_COLUMNS + FIT_COLUMNS, # -s3 with fit + 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 } # 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. diff --git a/ledsa/core/image_reading.py b/ledsa/core/image_reading.py index 7aa56d3..dd51536 100644 --- a/ledsa/core/image_reading.py +++ b/ledsa/core/image_reading.py @@ -38,6 +38,53 @@ def read_channel_data_from_img(filename: str, channel: int) -> np.ndarray: reader = _get_reader(filename, _read_channel_data_from_img_file, _read_channel_data_from_raw_file) return _read_or_raise(filename, reader, channel) + +def read_channel_data_and_mask_from_img(filename: str, channel: int) -> Tuple[np.ndarray, 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. + + 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. + + :param filename: The path of the image file to read. + :type filename: str + :param channel: The color channel to process. + :type channel: int + :return: The channel array and the boolean channel mask of the same shape. + :rtype: tuple[np.ndarray, np.ndarray] + :raises CorruptedImageError: If the image cannot be read or has an unsupported file extension. + """ + def read_img_file(img_filename, img_channel): + channel_array = _read_channel_data_from_img_file(img_filename, img_channel) + return channel_array, np.ones(channel_array.shape, dtype=bool) + + def read_raw_file(raw_filename, raw_channel): + img_array, filter_array = _read_img_array_from_raw_file(raw_filename, raw_channel) + channel_mask = get_bayer_channel_mask(filter_array, raw_channel) + return np.where(channel_mask, img_array, 0), channel_mask + + reader = _get_reader(filename, read_img_file, read_raw_file) + return _read_or_raise(filename, reader, channel) + + +def get_bayer_channel_mask(filter_array: np.ndarray, channel: int) -> np.ndarray: + """ + Returns a boolean mask of the Bayer filter positions that belong to the given color channel. + The Bayer array contains two green channels (1 and 3), one red (0) and one blue (2). + + :param filter_array: Bayer filter pattern of the raw image (rawpy raw_colors_visible). + :type filter_array: np.ndarray + :param channel: The color channel (0 red, 1 green, 2 blue). + :type channel: int + :return: Boolean mask, True where the pixel belongs to the channel. + :rtype: np.ndarray + """ + if channel == 1: + return (filter_array == 1) | (filter_array == 3) + return filter_array == channel + def read_img_array_from_img(filename: str, channel: int) -> np.ndarray: """ Returns a 2D array of the image for all color channels. diff --git a/ledsa/data_extraction/LEDAnalysisData.py b/ledsa/data_extraction/LEDAnalysisData.py index 6ecd04b..c2bcfca 100644 --- a/ledsa/data_extraction/LEDAnalysisData.py +++ b/ledsa/data_extraction/LEDAnalysisData.py @@ -21,6 +21,16 @@ class LEDAnalysisData: :ivar bgsub_sum_color_value: Integrated color value of the LED over the search area after subtraction of the local background estimated from the border of the search area. :vartype bgsub_sum_color_value: float + :ivar sidebg_sum_color_value: Integrated color value of the LED over the search area after subtraction of the + pedestal estimated in side bands next to the LED strip, at least a fraction of the plain sum. + :vartype sidebg_sum_color_value: float + :ivar sidebg_pedestal: Pedestal N * B subtracted for sidebg_sum_color_value (N channel pixels of the search + area, B mean per channel pixel in the side bands). + :vartype sidebg_pedestal: float + :ivar bgsub2_sum_color_value: Integrated color value of the LED over the search area minus N * b, with b the mean + over the channel pixels (including 0) of the border ring of the search area, at least a fraction of the + plain sum (sideband_floor). + :vartype bgsub2_sum_color_value: float :ivar fit_results: Fit results after fitting. :vartype fit_results: OptimizeResult :ivar fit_time: Time taken for fitting. @@ -45,6 +55,9 @@ def __init__(self, led_id, led_array, fit_leds): self.sum_color_value = None self.max_color_value = None self.bgsub_sum_color_value = None + self.sidebg_sum_color_value = None + self.sidebg_pedestal = None + self.bgsub2_sum_color_value = None self.fit_leds = fit_leds self.fit_results = None self.fit_time = None @@ -73,6 +86,8 @@ def get_main_data_string(self) -> str: out_str = f'{self.led_id:4d},{self.led_array:2d},' out_str += f'{self.sum_color_value:10.4e},{self.mean_color_value:10.4e},{self.max_color_value}' out_str += f',{self.bgsub_sum_color_value:10.4e}' + 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}' return out_str def get_fit_data_string(self) -> str: diff --git a/ledsa/data_extraction/sideband_background.py b/ledsa/data_extraction/sideband_background.py new file mode 100644 index 0000000..dcd7189 --- /dev/null +++ b/ledsa/data_extraction/sideband_background.py @@ -0,0 +1,186 @@ +""" +Pedestal (additive background) of the LED search areas, estimated in side bands next to the LED strips. + +The search area of an LED contains, besides the LED itself, an additive pedestal: stray light, light emitted or +scattered by the smoke and noise clipped at 0 after the black level subtraction of raw images. Under dense smoke the +LED dims while the pedestal remains, so the plain box sum (sum_col_val) overestimates the transmission. The border +ring of the search area (bgsub_sum_col_val) is not suited to measure the pedestal because it lies in the halo of the +LED itself. The side bands used here run parallel to the strip, at a distance of typically 25 to 60 px across the +strip and over +-R along the strip (R = search_area_radius), so they sample the scene next to the LED but not the LED. + +The strip direction is estimated locally for each LED from its nearest neighbours in the same LED array, so vertical, +horizontal and oblique strips are supported. Pixels close to any other search area are excluded from the bands, so +neighbouring strips do not enter the estimate. +""" +from typing import List, Optional, Sequence, Tuple + +import numpy as np + +DEFAULT_SIDEBAND_DISTANCES = (25, 60) +DEFAULT_SIDEBAND_FLOOR = 0.05 +NUM_DIRECTION_NEIGHBOURS = 4 + +# geometry of the side bands, cached per process because it only depends on the LED positions +_sideband_cache = {} + + +def estimate_strip_directions(search_areas: np.ndarray, led_arrays: Sequence[Sequence[int]], + num_neighbours: int = NUM_DIRECTION_NEIGHBOURS) -> np.ndarray: + """ + Estimate the local direction of the LED strip at each LED from its nearest neighbours in the same LED array. + + The direction is the principal axis of the positions of the LED and its num_neighbours nearest neighbours of the + same array. It is oriented such that the row component is positive (or, for horizontal strips, the column + component). + + :param search_areas: Array with rows (LED id, pixel row, pixel column), the row index equals the LED id. + :type search_areas: np.ndarray + :param led_arrays: LED ids of each LED array. + :type led_arrays: Sequence[Sequence[int]] + :param num_neighbours: Number of neighbours used to estimate the direction. + :type num_neighbours: int + :return: Array of shape (number of LEDs, 2) with the unit vector (row, column) along the strip, NaN for LEDs + that are not part of an array with at least two LEDs. + :rtype: np.ndarray + """ + directions = np.full((search_areas.shape[0], 2), np.nan) + for led_ids in led_arrays: + led_ids = np.atleast_1d(np.asarray(led_ids, dtype=int)) + if led_ids.size < 2: + continue + positions = search_areas[led_ids, 1:3].astype(np.float64) + num_points = min(num_neighbours, led_ids.size - 1) + 1 + for i, led_id in enumerate(led_ids): + distances = np.hypot(*(positions - positions[i]).T) + nearest = np.argsort(distances, kind='stable')[:num_points] + points = positions[nearest] - positions[nearest].mean(axis=0) + direction = np.linalg.svd(points, full_matrices=False)[2][0] + if direction[0] < 0 or (direction[0] == 0 and direction[1] < 0): + direction = -direction + directions[led_id] = direction + return directions + + +def build_sideband_pixels(search_areas: np.ndarray, directions: np.ndarray, img_shape: Tuple[int, int], + search_area_radius: int, distances: Tuple[float, float] = DEFAULT_SIDEBAND_DISTANCES, + exclusion_margin: Optional[int] = None) -> List[Optional[Tuple[np.ndarray, np.ndarray]]]: + """ + Determine the pixels of the two side bands of each LED. + + With t the coordinate along the strip and s the coordinate across the strip, both relative to the center of the + search area, a pixel belongs to the side bands if ``-R <= t < R`` and ``lo <= abs(s) < hi`` (with the same + half-open convention as the search area, i.e. ``lo <= s < hi`` or ``-hi <= s < -lo``). For a vertical strip this is + rows r-R..r+R-1 and columns c+lo..c+hi-1 and c-hi..c-lo-1. Pixels outside the image and pixels within + R + exclusion_margin (Chebyshev distance) of the center of any other search area are excluded. + + :param search_areas: Array with rows (LED id, pixel row, pixel column), the row index equals the LED id. + :type search_areas: np.ndarray + :param directions: Unit vectors along the strip per LED, see estimate_strip_directions. + :type directions: np.ndarray + :param img_shape: Shape of the image. + :type img_shape: tuple[int, int] + :param search_area_radius: Radius R of the search areas. + :type search_area_radius: int + :param distances: Inner and outer distance (lo, hi) of the side bands across the strip in pixels. + :type distances: tuple[float, float] + :param exclusion_margin: Margin added to R for the exclusion of other search areas. Defaults to R. + :type exclusion_margin: int or None + :return: Per LED id a tuple of row and column indices of the side band pixels, None if the strip direction of + the LED is unknown. + :rtype: list[tuple[np.ndarray, np.ndarray] or None] + """ + radius = int(search_area_radius) + inner, outer = float(distances[0]), float(distances[1]) + if exclusion_margin is None: + exclusion_margin = radius + exclusion = radius + int(exclusion_margin) + reach = int(np.ceil(np.hypot(outer, radius))) + 1 + d_row, d_col = np.mgrid[-reach:reach + 1, -reach:reach + 1] + centers = search_areas[:, 1:3].astype(np.int64) + sideband_pixels = [None] * search_areas.shape[0] + for led_id in range(search_areas.shape[0]): + along = directions[led_id] + if not np.all(np.isfinite(along)): + continue + across = np.array([-along[1], along[0]]) + # rounding avoids that pixels exactly on a band border flip due to floating point errors + t = np.round(d_row * along[0] + d_col * along[1], 6) + s = np.round(d_row * across[0] + d_col * across[1], 6) + in_band = (t >= -radius) & (t < radius) & (((s >= inner) & (s < outer)) | ((s >= -outer) & (s < -inner))) + rows = centers[led_id, 0] + d_row[in_band] + cols = centers[led_id, 1] + d_col[in_band] + keep = (rows >= 0) & (rows < img_shape[0]) & (cols >= 0) & (cols < img_shape[1]) + nearby = np.abs(centers - centers[led_id]).max(axis=1) <= reach + exclusion + nearby[led_id] = False + for row_other, col_other in centers[nearby]: + keep &= ~((rows >= row_other - exclusion) & (rows < row_other + exclusion) & + (cols >= col_other - exclusion) & (cols < col_other + exclusion)) + sideband_pixels[led_id] = (rows[keep].astype(np.intp), cols[keep].astype(np.intp)) + return sideband_pixels + + +def get_sideband_pixels(search_areas: np.ndarray, led_arrays: Sequence[Sequence[int]], img_shape: Tuple[int, int], + search_area_radius: int, distances: Tuple[float, float] = DEFAULT_SIDEBAND_DISTANCES, + exclusion_margin: Optional[int] = None) -> List[Optional[Tuple[np.ndarray, np.ndarray]]]: + """ + Cached combination of estimate_strip_directions and build_sideband_pixels. The geometry only depends on the + LED positions and is therefore computed once per process and reused for all images. + + :return: See build_sideband_pixels. + """ + key = (search_areas.tobytes(), search_areas.shape, + tuple(tuple(int(i) for i in np.atleast_1d(led_ids)) for led_ids in led_arrays), + tuple(img_shape), int(search_area_radius), tuple(float(d) for d in distances), exclusion_margin) + if key not in _sideband_cache: + _sideband_cache.clear() + directions = estimate_strip_directions(search_areas, led_arrays) + _sideband_cache[key] = build_sideband_pixels(search_areas, directions, img_shape, search_area_radius, + distances, exclusion_margin) + return _sideband_cache[key] + + +def estimate_sideband_background(data: np.ndarray, channel_mask: Optional[np.ndarray], + sideband_pixels: Optional[Tuple[np.ndarray, np.ndarray]]) -> float: + """ + Estimate the pedestal per channel pixel as the mean over the channel pixels in the side bands. + + Pixels of the channel with value 0 are included: the black level subtraction clips noise at 0 and the box sum + contains these clipped pixels in the same way, so the mean is the matching estimator for the pedestal in the box. + Foreign Bayer channels are excluded via the channel mask, not via the pixel value. + + :param data: Channel array of the image. + :type data: np.ndarray + :param channel_mask: Boolean mask of the channel pixels, None if all pixels belong to the channel. + :type channel_mask: np.ndarray or None + :param sideband_pixels: Row and column indices of the side band pixels of the LED. + :type sideband_pixels: tuple[np.ndarray, np.ndarray] or None + :return: Mean pedestal per channel pixel, NaN if the side bands contain no channel pixel. + :rtype: float + """ + if sideband_pixels is None or sideband_pixels[0].size == 0: + return np.nan + values = data[sideband_pixels] + if channel_mask is not None: + values = values[channel_mask[sideband_pixels]] + if values.size == 0: + return np.nan + return float(np.mean(values, dtype=np.float64)) + + +def subtract_pedestal(sum_color_value: float, pedestal: float, floor: float = DEFAULT_SIDEBAND_FLOOR) -> float: + """ + Subtract the pedestal N * B from the box sum S with a lower limit: max(S - N * B, floor * S). + + Without the lower limit the optical depth would be unbounded where the LED signal vanishes in the pedestal; + values where the limit applies carry no information on the extinction and should be treated as invalid. + + :param sum_color_value: Box sum S. + :type sum_color_value: float + :param pedestal: Pedestal N * B of the box. + :type pedestal: float + :param floor: Lower limit as fraction of S. + :type floor: float + :return: Pedestal corrected box sum, NaN if the pedestal is NaN. + :rtype: float + """ + return float(np.maximum(sum_color_value - pedestal, floor * sum_color_value)) diff --git a/ledsa/data_extraction/step_3_functions.py b/ledsa/data_extraction/step_3_functions.py index 2b30e72..ab8edd4 100644 --- a/ledsa/data_extraction/step_3_functions.py +++ b/ledsa/data_extraction/step_3_functions.py @@ -1,7 +1,7 @@ import os import re import time -from typing import List, Tuple +from typing import List, Optional, Tuple import numpy as np import scipy.optimize @@ -9,9 +9,12 @@ 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_from_img +from ledsa.core.image_reading import read_channel_data_and_mask_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_SIDEBAND_FLOOR, + estimate_sideband_background, get_sideband_pixels, + subtract_pedestal) def generate_analysis_data(img_filename: str, channel: int, search_areas: np.ndarray, line_indices: List[List[int]], @@ -39,7 +42,7 @@ 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 = read_channel_data_from_img(file_path, channel=channel) + data, channel_mask = read_channel_data_and_mask_from_img(file_path, channel=channel) search_area_radius = int(conf['find_search_areas']['search_area_radius']) img_analysis_data = [] @@ -48,16 +51,70 @@ def generate_analysis_data(img_filename: str, channel: int, search_areas: np.nda search_area_radius, fit_leds) return analysis_res + sideband_distances, sideband_floor, sideband_exclusion_margin = get_sideband_parameters(conf) + sideband_pixels = get_sideband_pixels(search_areas, _get_unmerged_led_arrays(conf, line_indices), data.shape, + search_area_radius, sideband_distances, sideband_exclusion_margin) + num_of_arrays = len(line_indices) for led_array_idx in range(num_of_arrays): for iled in line_indices[led_array_idx]: if iled % (int(conf['analyse_photo']['num_skip_leds']) + 1) == 0: led_analysis_data = _generate_led_analysis_data(conf, channel, data, debug, iled, img_filename, - led_array_idx, search_areas, search_area_radius, fit_leds) + led_array_idx, search_areas, search_area_radius, fit_leds, + channel_mask=channel_mask, + sideband_pixels=sideband_pixels[iled], + sideband_floor=sideband_floor) img_analysis_data.append(led_analysis_data) return img_analysis_data +def get_sideband_parameters(conf: ConfigData) -> Tuple[Tuple[float, float], float, Optional[int]]: + """ + Read the parameters of the side band pedestal estimation from the config. Config files written before the + side band estimation was added do not contain the options, the defaults are used then. + + :param conf: Configuration data. + :type conf: ConfigData + :return: Distances (inner, outer) of the side bands across the strip in pixels, lower limit of the corrected + sum as fraction of the box sum and margin for the exclusion of other search areas (None: R). + :rtype: tuple[tuple[float, float], float, int or None] + """ + distances = DEFAULT_SIDEBAND_DISTANCES + floor = DEFAULT_SIDEBAND_FLOOR + exclusion_margin = None + if conf is not None: + if conf.get('analyse_photo', 'sideband_distances') != 'None': + distances = tuple(float(d) for d in conf.get('analyse_photo', 'sideband_distances').split()) + if len(distances) != 2 or not 0 <= distances[0] < distances[1]: + raise ValueError(f'sideband_distances must be two distances 0 <= inner < outer, got {distances}.') + if conf.get('analyse_photo', 'sideband_floor') != 'None': + floor = float(conf.get('analyse_photo', 'sideband_floor')) + if conf.get('analyse_photo', 'sideband_exclusion_margin') != 'None': + exclusion_margin = int(conf.get('analyse_photo', 'sideband_exclusion_margin')) + return distances, floor, exclusion_margin + + +def _get_unmerged_led_arrays(conf: ConfigData, line_indices: List[List[int]]) -> List[List[int]]: + """ + Return the LED arrays as found in step 2. If LED arrays are merged for the analysis, the merged arrays may + contain several strips, which would spoil the estimate of the local strip direction. + + :param conf: Configuration data. + :type conf: ConfigData + :param line_indices: LED ids of the (possibly merged) LED arrays used in the analysis. + :type line_indices: List[List[int]] + :return: LED ids of the LED arrays as found in step 2. + :rtype: List[List[int]] + """ + if conf is None or conf['analyse_positions']['merge_led_array_indices'] == 'None': + return line_indices + led_arrays = [] + for i in range(int(conf['analyse_positions']['num_arrays'])): + file_path = os.path.join('analysis', f'led_array_indices_{i:03}.csv') + led_arrays.append(read_table(file_path, dtype='int', silent=True)) + return led_arrays + + def create_fit_result_file(img_data: List[LEDAnalysisData], img_id: int, channel: int) -> None: # TODO: rename because misleading """ Create a result file for a single image, containing the pixel values and, if applicable, the fit results of all LEDs. @@ -107,7 +164,10 @@ def find_and_save_not_analysed_imgs(channel: int) -> None: _save_list_of_remaining_imgs_needed_to_be_processed(remaining_imgs) -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) -> LEDAnalysisData: +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, + sideband_pixels: Optional[Tuple[np.ndarray, np.ndarray]] = None, + sideband_floor: float = DEFAULT_SIDEBAND_FLOOR) -> LEDAnalysisData: """ Generate analysis data for a specific LED. @@ -131,6 +191,13 @@ 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 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 + :param sideband_pixels: Row and column indices of the side band pixels of the LED, None if not available. + :type sideband_pixels: tuple[np.ndarray, np.ndarray] or None + :param sideband_floor: Lower limit of the side band and ring corrected sums as fraction of the box sum. + :type sideband_floor: float :return: Analysis data for the LED. :rtype: LEDAnalysisData """ @@ -162,6 +229,22 @@ 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))) + # 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)) + + # border ring as for bgsub, but with the estimator of the side bands: mean over the channel pixels including 0; + # same lower limit as for the side band correction + ring_background = _estimate_local_background_mean(data[search_area], box_mask) + led_data.bgsub2_sum_color_value = subtract_pedestal(float(led_data.sum_color_value), + num_channel_pixels * ring_background, sideband_floor) + + # b = mean over the channel pixels in the side bands next to the LED strip + sideband_background = estimate_sideband_background(data, channel_mask, sideband_pixels) + led_data.sidebg_pedestal = num_channel_pixels * sideband_background + led_data.sidebg_sum_color_value = subtract_pedestal(float(led_data.sum_color_value), led_data.sidebg_pedestal, + sideband_floor) + return led_data @@ -192,6 +275,38 @@ def _estimate_local_background(search_area_data: np.ndarray, border_width: int = return float(np.median(border_pixels)) +def _estimate_local_background_mean(search_area_data: np.ndarray, channel_mask: Optional[np.ndarray] = None, + border_width: int = 2) -> float: + """ + Estimate the local background of a search area as the mean over the channel pixels of its border ring. + + Unlike _estimate_local_background, pixels with value 0 are kept if they belong to the color channel: they are + noise clipped at 0 after the black level subtraction and are part of the box sum in the same way. Pixels of + foreign Bayer channels are excluded via the channel mask. The result is meant to be subtracted as N * b from the + box sum (N channel pixels in the box), which is linear in the pixel values, instead of pixel-wise with clipping. + Note that the border ring lies in the halo of the LED for small search areas. + + :param search_area_data: Part of the image where the LED is located. + :type search_area_data: np.ndarray + :param channel_mask: Boolean mask of the channel pixels in the search area, None if all pixels belong to the + channel. + :type channel_mask: np.ndarray or None + :param border_width: Width of the border ring in pixels. + :type border_width: int + :return: Mean background per channel pixel, NaN if the border ring contains no channel pixel. + :rtype: float + """ + if search_area_data.shape[0] <= 2 * border_width or search_area_data.shape[1] <= 2 * border_width: + return np.nan + border_mask = np.ones(search_area_data.shape, dtype=bool) + border_mask[border_width:-border_width, border_width:-border_width] = False + if channel_mask is not None: + border_mask &= channel_mask + if not border_mask.any(): + return np.nan + return float(np.mean(search_area_data[border_mask], dtype=np.float64)) + + def _save_results_in_file(channel: int, img_data: LEDAnalysisData, img_filename: str, img_id: str, img_infos: np.ndarray, basename: str) -> None: """ Save analysis results to a file. @@ -333,6 +448,7 @@ def _create_header(channel: int, img_id: str, img_filename: str, img_infos: np.n out_str += f"channel = {channel}, " out_str += f"time[s] = {img_infos[int(img_id) - 1][3]}\n" out_str += "# id,line,sum_col_value,average_col_value,max_col_value,bgsub_sum_col_value" + out_str += ",sidebg_sum_col_value,sidebg_pedestal,bgsub2_sum_col_value" if fit_leds: out_str += ",led_center_x, led_center_y" out_str += ",x,y,dx,dy,A,alpha,wx,wy,fit_success,fit_fun,fit_nfev,fit_time" diff --git a/ledsa/tests/AcceptanceTests/10_test_sideband_background.robot b/ledsa/tests/AcceptanceTests/10_test_sideband_background.robot new file mode 100644 index 0000000..58f7fac --- /dev/null +++ b/ledsa/tests/AcceptanceTests/10_test_sideband_background.robot @@ -0,0 +1,89 @@ +*** Settings *** +Resource global_keywords.resource + +Force Tags analysis background_subtraction sideband + +*** Variables *** +${BINARY} analysis${/}channel0${/}all_parameters.h5 +${EXTINCTION_DIR} analysis${/}extinction_coefficients${/}linear +${SUM_RESULT} ${EXTINCTION_DIR}${/}extinction_coefficients_linear_channel_0_sum_col_val_led_array_0.csv +${SUM_REFERENCE} ${EXTINCTION_DIR}${/}extinction_coefficients_linear_channel_0_sum_col_val_led_array_0_sidebg_ref.csv +${BGSUB_RESULT} ${EXTINCTION_DIR}${/}extinction_coefficients_linear_channel_0_bgsub_sum_col_val_led_array_0.csv +${BGSUB_REFERENCE} ${EXTINCTION_DIR}${/}extinction_coefficients_linear_channel_0_bgsub_sum_col_val_led_array_0_sidebg_ref.csv + +*** Test Cases *** +Default Results Are Not Changed By The Sideband Columns + [Documentation] Results with the existing reference properties sum_col_val and bgsub_sum_col_val must be + ... identical whether the step 3 files contain the side band columns or not (data format + ... before the side band pedestal was added). The binary is removed before each comparison + ... run, because it is only created from the step 3 files if it does not exist. + Change Directory ${WORKDIR} + Execute Ledsa -s3_fast + Remove File ${BINARY} + Create And Fill Config Analysis linear + Execute Ledsa --analysis + Copy File ${SUM_RESULT} ${SUM_REFERENCE} + Create And Fill Config Analysis linear bgsub_sum_col_val + Execute Ledsa --analysis + Copy File ${BGSUB_RESULT} ${BGSUB_REFERENCE} + Strip Bgsub Column From Led Position Files num_columns=6 + Remove File ${BINARY} + Execute Ledsa --analysis + ${identical} = Check Extinction Coefficient Files Are Identical ${BGSUB_RESULT} ${BGSUB_REFERENCE} + Should Be True ${identical} + Create And Fill Config Analysis linear + Execute Ledsa --analysis + ${identical} = Check Extinction Coefficient Files Are Identical ${SUM_RESULT} ${SUM_REFERENCE} + Should Be True ${identical} + +Sidebg Reference Property Matches Ground Truth On Clean Images + [Documentation] Without background the side bands contain nothing to subtract and the sidebg property + ... reproduces the ground truth like the default property. + Run Sideband Analysis sideband_clean 0 + Check Sidebg Results 2 + Check Sidebg Results 3 + Check Sidebg Results 4 + +Background Light Biases Sum But Not Sidebg + [Documentation] With a constant additive background the default property underestimates the extinction + ... while the sidebg property stays close to the ground truth. + Run Sideband Analysis sideband_background 15 + Check Background Results 2 + Check Background Results 3 + Check Background Results 4 + +*** Keywords *** +Run Sideband Analysis + [Arguments] ${dir} ${background} + [Documentation] The test images are only 50 px wide, so the side bands are placed 12 to 20 px across + ... the strip instead of the default 25 to 60 px. + Create Directory ${WORKDIR}${/}${dir} + Change Directory ${WORKDIR}${/}${dir} + Create Test Data background=${background} + Create Config + Set Config Option analyse_photo sideband_distances 12 20 + Execute Ledsa -s1 + Execute Ledsa -s2 + Execute Ledsa --coordinates + Execute Ledsa -s3_fast + Create And Fill Config Analysis linear + Execute Ledsa --analysis + Create And Fill Config Analysis linear sidebg_sum_col_val + Execute Ledsa --analysis + +Check Sidebg Results + [Arguments] ${image_id} + ${rmse} = Check Input Vs Computed Extinction Coefficients ${image_id} linear + ... reference_property=sidebg_sum_col_val + Should Be True ${rmse} < 0.05 + +Check Background Results + [Arguments] ${image_id} + ${rmse_sum} = Check Input Vs Computed Extinction Coefficients ${image_id} linear + ${rmse_sidebg} = Check Input Vs Computed Extinction Coefficients ${image_id} linear + ... reference_property=sidebg_sum_col_val + # same slack as for bgsub in suite 08: the test data are random and the background clips the LED peaks of + # the 8 bit images at 255; in addition the side band correction has to remove most of the bias of the sum + Should Be True ${rmse_sum} > 0.05 + Should Be True ${rmse_sidebg} < 0.06 + Should Be True ${rmse_sidebg} < 0.5 * ${rmse_sum} diff --git a/ledsa/tests/UnitTests/test_background_subtraction.py b/ledsa/tests/UnitTests/test_background_subtraction.py index 829d805..7363909 100644 --- a/ledsa/tests/UnitTests/test_background_subtraction.py +++ b/ledsa/tests/UnitTests/test_background_subtraction.py @@ -90,16 +90,31 @@ def test_sum_col_value_is_unchanged_by_new_quantity(self): def test_csv_serialization_contains_bgsub_value(self): led_data = extract_led_data(make_synthetic_image(background=10.0)) main_data_fields = led_data.get_main_data_string().split(',') - assert len(main_data_fields) == 6 + assert len(main_data_fields) >= 6 # further background columns are appended after bgsub assert float(main_data_fields[5]) == pytest.approx(led_data.bgsub_sum_color_value, rel=1e-4) class TestColumnNames: - def _write_led_positions_csv(self, tmp_path, num_columns): + def _write_led_positions_csv(self, tmp_path, num_columns, header='# header'): channel_dir = tmp_path / 'analysis' / 'channel0' channel_dir.mkdir(parents=True) row = ','.join(str(float(i)) for i in range(num_columns)) - (channel_dir / '1_led_positions.csv').write_text(f'# header\n{row}\n{row}\n', encoding='utf-8') + (channel_dir / '1_led_positions.csv').write_text(f'{header}\n{row}\n{row}\n', encoding='utf-8') + + def test_csv_written_by_step_3_without_fit(self, tmp_path, monkeypatch): + from ledsa.core.file_handling import _get_column_names + from ledsa.data_extraction.step_3_functions import _create_header + img_infos = np.array([['1', 'img.jpg', '12:00:00', '0']]) + header = _create_header(0, '1', 'img.jpg', img_infos, 'root', False).strip() + num_columns = len(header.split('\n')[1].split(',')) + self._write_led_positions_csv(tmp_path, num_columns, header) + monkeypatch.chdir(tmp_path) + columns = _get_column_names(0) + assert columns[:7] == ["img_id", "led_id", "led_array_id", "sum_col_val", "mean_col_val", + "max_col_val", "bgsub_sum_col_val"] + assert "sidebg_sum_col_val" in columns and "sidebg_pedestal" in columns + assert columns[-2:] == ["width", "height"] + assert len(columns) == num_columns + 3 def test_new_csv_with_bgsub_column(self, tmp_path, monkeypatch): from ledsa.core.file_handling import _get_column_names @@ -125,6 +140,8 @@ def test_fitted_csv_written_by_step_3(self, tmp_path, monkeypatch): led_data = LEDAnalysisData(3, 0, True) led_data.sum_color_value, led_data.mean_color_value, led_data.max_color_value = 41210.0, 103.0, 232 led_data.bgsub_sum_color_value = 40000.0 + led_data.sidebg_sum_color_value, led_data.sidebg_pedestal = 39000.0, 2210.0 + led_data.bgsub2_sum_color_value = 39500.0 led_data.led_center_x, led_data.led_center_y = 50.521, 60.25 led_data.fit_results = SimpleNamespace(x=np.arange(1.5, 9.5), success=True, fun=0.5, nfev=120) led_data.fit_time = 0.25 @@ -138,11 +155,26 @@ def test_fitted_csv_written_by_step_3(self, tmp_path, monkeypatch): named = dict(zip(columns[1:], values)) assert named['max_col_val'] == 232 assert named['bgsub_sum_col_val'] == pytest.approx(40000.0) + assert named['sidebg_pedestal'] == pytest.approx(2210.0) + assert named['bgsub2_sum_col_val'] == pytest.approx(39500.0) assert named['led_center_x'] == pytest.approx(50.521) assert named['A'] == pytest.approx(5.5) assert named['fit_nfev'] == 120 assert named['fit_time'] == pytest.approx(0.25) + @pytest.mark.parametrize('num_columns, last_main', [(6, 'bgsub_sum_col_val'), (9, 'bgsub2_sum_col_val'), + (20, 'bgsub_sum_col_val'), (23, 'bgsub2_sum_col_val')]) + def test_layouts(self, tmp_path, monkeypatch, num_columns, last_main): + from ledsa.core.file_handling import _get_column_names + self._write_led_positions_csv(tmp_path, num_columns) + monkeypatch.chdir(tmp_path) + columns = _get_column_names(0) + assert len(columns) == num_columns + 1 + 2 + main = columns[:columns.index('led_center_x')] if 'led_center_x' in columns else columns[:-2] + assert main[-1] == last_main + if num_columns > 9: + assert columns[-3:] == ['fit_time', 'width', 'height'] + def test_fitted_csv_without_delimiter_raises_with_hint(self, tmp_path, monkeypatch): # fitted files of this branch written before the delimiter fix: bgsub_sum_col_value and led_center_x glued from ledsa.core.file_handling import _get_column_names @@ -177,7 +209,7 @@ def test_oldest_layouts(self, tmp_path, monkeypatch, num_columns, first_names): assert columns[:4] == first_names assert len(columns) == num_columns + 1 + 2 - @pytest.mark.parametrize('num_columns', [7, 19, 21]) + @pytest.mark.parametrize('num_columns', [7, 19, 21, 22]) def test_unknown_layout_raises(self, tmp_path, monkeypatch, num_columns): from ledsa.core.file_handling import _get_column_names self._write_led_positions_csv(tmp_path, num_columns) diff --git a/ledsa/tests/UnitTests/test_ring_background_mean.py b/ledsa/tests/UnitTests/test_ring_background_mean.py new file mode 100644 index 0000000..4adcb19 --- /dev/null +++ b/ledsa/tests/UnitTests/test_ring_background_mean.py @@ -0,0 +1,102 @@ +import numpy as np +import pytest + +from ledsa.data_extraction.step_3_functions import (_estimate_local_background, _estimate_local_background_mean, + _generate_led_analysis_data) + +SEARCH_AREA_RADIUS = 8 +CENTER = (40, 40) + + +def make_image(transmission=1.0, pedestal=0.0, noise_sigma=0.0, clip=False, bayer_mask=False, seed=0, + amplitude=4000.0, sigma=1.5): + """Single Gaussian LED plus constant pedestal and noise, optionally clipped at 0 and with an RGGB green mask.""" + rng = np.random.default_rng(seed) + rows, cols = np.mgrid[0:80, 0:80] + led = amplitude * transmission * np.exp(-((rows - CENTER[0]) ** 2 + (cols - CENTER[1]) ** 2) / (2 * sigma ** 2)) + img = led + pedestal + rng.normal(0.0, noise_sigma, led.shape) + if clip: + img = np.clip(img, 0.0, None) + mask = (rows + cols) % 2 == 1 if bayer_mask else np.ones(img.shape, dtype=bool) + return np.where(mask, img, 0.0), mask + + +def analyse(img, mask): + search_areas = np.array([[0, CENTER[0], CENTER[1]]]) + return _generate_led_analysis_data(None, 1, img, False, 0, 'test_img', 0, search_areas, SEARCH_AREA_RADIUS, + fit_leds=False, channel_mask=mask) + + +def box(array): + r = SEARCH_AREA_RADIUS + return array[CENTER[0] - r:CENTER[0] + r, CENTER[1] - r:CENTER[1] + r] + + +class TestRingBackgroundMean: + def test_constant_image(self): + assert _estimate_local_background_mean(np.full((16, 16), 7.5)) == pytest.approx(7.5) + + def test_foreign_channel_pixels_are_excluded_by_mask(self): + img, mask = make_image(pedestal=7.5, bayer_mask=True, amplitude=0.0) + assert _estimate_local_background_mean(box(img), box(mask)) == pytest.approx(7.5) + + def test_channel_pixels_with_value_zero_are_included(self): + data = np.full((16, 16), 4.0) + data[::2, :] = 0.0 # channel pixels clipped at 0 + ring = np.ones((16, 16), dtype=bool) + ring[2:-2, 2:-2] = False + assert _estimate_local_background_mean(data) == pytest.approx(data[ring].mean()) + assert _estimate_local_background(data) == pytest.approx(4.0) # old estimator ignores the zeros + + def test_too_small_or_empty_gives_nan(self): + assert np.isnan(_estimate_local_background_mean(np.ones((4, 4)))) + assert np.isnan(_estimate_local_background_mean(np.ones((16, 16)), np.zeros((16, 16), dtype=bool))) + + +class TestBgsub2Sum: + @pytest.mark.parametrize('bayer_mask', [False, True]) + def test_subtraction_is_linear(self, bayer_mask): + img, mask = make_image(pedestal=3.0, noise_sigma=2.0, clip=True, bayer_mask=bayer_mask) + led_data = analyse(img, mask) + num_channel_pixels = np.count_nonzero(box(mask)) + expected = box(img).sum() - num_channel_pixels * _estimate_local_background_mean(box(img), box(mask)) + assert led_data.bgsub2_sum_color_value == pytest.approx(expected) + + def test_floor_limits_the_corrected_sum(self): + # LED shifted to the corner of the box: its light falls into the ring and N * b exceeds the box sum + img, mask = make_image(amplitude=0.0) + row, col = CENTER[0] - SEARCH_AREA_RADIUS + 1, CENTER[1] - SEARCH_AREA_RADIUS + 1 + img[row - 1:row + 2, col - 1:col + 2] = 1000.0 + led_data = analyse(img, mask) + num_channel_pixels = np.count_nonzero(box(mask)) + assert box(img).sum() - num_channel_pixels * _estimate_local_background_mean(box(img), box(mask)) < 0 + assert led_data.bgsub2_sum_color_value == pytest.approx(0.05 * box(img).sum()) + + def test_existing_values_are_unchanged(self): + img, mask = make_image(pedestal=3.0, noise_sigma=2.0, clip=True, bayer_mask=True) + led_data = analyse(img, mask) + reference = _generate_led_analysis_data(None, 1, img, False, 0, 'test_img', 0, np.array([[0, *CENTER]]), + SEARCH_AREA_RADIUS, fit_leds=False) + assert led_data.sum_color_value == reference.sum_color_value + assert led_data.bgsub_sum_color_value == reference.bgsub_sum_color_value + fields = led_data.get_main_data_string().split(',') + assert ','.join(fields[:6]) == ','.join(reference.get_main_data_string().split(',')[:6]) + assert float(fields[8]) == pytest.approx(led_data.bgsub2_sum_color_value, rel=1e-4) + + def test_optical_depth_bias_with_clipped_noise(self): + # LED without halo in the ring: the mean estimator removes the pedestal including the clipped noise, + # the median of the pixels > 0 with pixel-wise clipping (bgsub) leaves a bias + transmission = 0.01 + args = dict(pedestal=2.0, noise_sigma=2.0, clip=True, bayer_mask=True) + true_depth = -np.log(transmission) + depth = {'sum': [], 'bgsub': [], 'bgsub2': []} + for seed in range(20): + ref = analyse(*make_image(transmission=1.0, seed=2 * seed, **args)) + img = analyse(*make_image(transmission=transmission, seed=2 * seed + 1, **args)) + depth['sum'].append(-np.log(img.sum_color_value / ref.sum_color_value)) + depth['bgsub'].append(-np.log(img.bgsub_sum_color_value / ref.bgsub_sum_color_value)) + depth['bgsub2'].append(-np.log(img.bgsub2_sum_color_value / ref.bgsub2_sum_color_value)) + bias = {key: np.mean(values) / true_depth - 1 for key, values in depth.items()} + assert abs(bias['bgsub2']) < 0.02 + assert bias['bgsub'] < -0.02 + assert bias['sum'] < -0.1 diff --git a/ledsa/tests/UnitTests/test_sideband_background.py b/ledsa/tests/UnitTests/test_sideband_background.py new file mode 100644 index 0000000..820c5c5 --- /dev/null +++ b/ledsa/tests/UnitTests/test_sideband_background.py @@ -0,0 +1,307 @@ +import numpy as np +import pytest + +from ledsa.data_extraction.sideband_background import (build_sideband_pixels, estimate_sideband_background, + estimate_strip_directions, get_sideband_pixels, + subtract_pedestal) +import ledsa.core.image_reading as image_reading +from ledsa.core.image_reading import get_bayer_channel_mask +from ledsa.data_extraction.step_3_functions import (_generate_led_analysis_data, _get_unmerged_led_arrays, + get_sideband_parameters) + +SEARCH_AREA_RADIUS = 8 +LED_SPACING = 20 +LED_AMPLITUDE = 400.0 +LED_SIGMA = 1.5 +IMG_SIZE = 240 +SIDEBAND_DISTANCES = (25, 60) + + +def make_strip(angle_deg=0.0, num_leds=7, center=(120, 120)): + """LED positions (row, col) along a straight strip; angle 0 is a vertical strip (along the rows).""" + along = np.array([np.cos(np.radians(angle_deg)), np.sin(np.radians(angle_deg))]) + offsets = (np.arange(num_leds) - num_leds // 2) * LED_SPACING + positions = np.rint(np.array(center) + offsets[:, None] * along).astype(int) + return np.column_stack([np.arange(num_leds), positions]) + + +def make_image(search_areas, transmission=1.0, pedestal=0.0, noise_sigma=0.0, bayer_mask=False, clip=False, + seed=0, amplitude=LED_AMPLITUDE): + """Gaussian LEDs scaled by the transmission plus a constant pedestal and Gaussian noise. + + With clip=True the image is clipped at 0 like LEDSA does after the black level subtraction of raw files. + With bayer_mask=True only the green pixels of an RGGB pattern belong to the channel, all other pixels are 0. + Returns the channel array, the channel mask and the pure LED image (without pedestal and noise). + """ + rng = np.random.default_rng(seed) + rows, cols = np.mgrid[0:IMG_SIZE, 0:IMG_SIZE] + led_img = np.zeros((IMG_SIZE, IMG_SIZE)) + for _, row, col in search_areas: + led_img += amplitude * transmission * np.exp(-((rows - row) ** 2 + (cols - col) ** 2) / (2 * LED_SIGMA ** 2)) + img = led_img + pedestal + rng.normal(0.0, noise_sigma, led_img.shape) + if clip: + img = np.clip(img, 0.0, None) + if bayer_mask: + mask = (rows + cols) % 2 == 1 + else: + mask = np.ones(img.shape, dtype=bool) + img = np.where(mask, img, 0.0) + led_img = np.where(mask, led_img, 0.0) + return img, mask, led_img + + +def box_sum(img, row, col, radius=SEARCH_AREA_RADIUS): + return float(np.sum(img[row - radius:row + radius, col - radius:col + radius])) + + +def led_data_for(img, mask, search_areas, led_id, floor=0.05): + sideband_pixels = get_sideband_pixels(search_areas, [search_areas[:, 0]], img.shape, SEARCH_AREA_RADIUS, + SIDEBAND_DISTANCES) + return _generate_led_analysis_data(None, 1, img, False, led_id, 'test_img', 0, search_areas, SEARCH_AREA_RADIUS, + fit_leds=False, channel_mask=mask, sideband_pixels=sideband_pixels[led_id], + sideband_floor=floor) + + +class TestStripDirection: + @pytest.mark.parametrize('angle', [0.0, 90.0, 30.0, -45.0]) + def test_direction_follows_strip(self, angle): + search_areas = make_strip(angle) + directions = estimate_strip_directions(search_areas, [search_areas[:, 0]]) + expected = np.array([np.cos(np.radians(angle)), np.sin(np.radians(angle))]) + for direction in directions: + # direction is defined up to the sign; rounding of the LED positions tilts it slightly + assert abs(np.dot(direction, expected)) == pytest.approx(1.0, abs=2e-3) + + def test_led_without_array_has_no_direction(self): + search_areas = make_strip(0.0) + directions = estimate_strip_directions(search_areas, [search_areas[:3, 0]]) + assert np.all(np.isfinite(directions[:3])) + assert np.all(np.isnan(directions[3:])) + + +class TestSidebandPixels: + def test_vertical_strip_matches_axis_aligned_bands(self): + # same pixels as the axis-aligned bands of the H2 evaluation: rows r-R..r+R-1, columns c+25..c+59 and + # c-60..c-26 + search_areas = make_strip(0.0) + directions = np.tile([1.0, 0.0], (len(search_areas), 1)) + pixels = build_sideband_pixels(search_areas, directions, (IMG_SIZE, IMG_SIZE), SEARCH_AREA_RADIUS, + SIDEBAND_DISTANCES) + led_id, row, col = search_areas[3] + expected = np.zeros((IMG_SIZE, IMG_SIZE), dtype=bool) + expected[row - SEARCH_AREA_RADIUS:row + SEARCH_AREA_RADIUS, col + 25:col + 60] = True + expected[row - SEARCH_AREA_RADIUS:row + SEARCH_AREA_RADIUS, col - 60:col - 25] = True + selected = np.zeros_like(expected) + selected[pixels[led_id]] = True + np.testing.assert_array_equal(selected, expected) + + @pytest.mark.parametrize('angle', [0.0, 90.0, 30.0]) + def test_bands_do_not_touch_search_areas(self, angle): + search_areas = make_strip(angle) + pixels = get_sideband_pixels(search_areas, [search_areas[:, 0]], (IMG_SIZE, IMG_SIZE), SEARCH_AREA_RADIUS, + SIDEBAND_DISTANCES) + for _, row, col in search_areas: + box = np.zeros((IMG_SIZE, IMG_SIZE), dtype=bool) + box[row - SEARCH_AREA_RADIUS:row + SEARCH_AREA_RADIUS, col - SEARCH_AREA_RADIUS:col + SEARCH_AREA_RADIUS] = True + for led_pixels in pixels: + assert led_pixels[0].size > 0 + assert not box[led_pixels].any() + + def test_neighbouring_strip_is_excluded(self): + strip = make_strip(0.0, center=(120, 80)) + neighbour = make_strip(0.0, center=(120, 120)) # 40 px to the right, inside the right side band + neighbour[:, 0] += len(strip) + search_areas = np.vstack([strip, neighbour]) + exclusion = 2 * SEARCH_AREA_RADIUS + pixels = get_sideband_pixels(search_areas, [strip[:, 0], neighbour[:, 0]], (IMG_SIZE, IMG_SIZE), + SEARCH_AREA_RADIUS, SIDEBAND_DISTANCES) + rows, cols = pixels[3] + assert np.all((cols < 120 - exclusion) | (cols >= 120 + exclusion)) + assert np.any(cols >= 120 + exclusion) + + def test_bands_are_clipped_at_the_image_border(self): + search_areas = make_strip(0.0, center=(120, 30)) + pixels = get_sideband_pixels(search_areas, [search_areas[:, 0]], (IMG_SIZE, IMG_SIZE), SEARCH_AREA_RADIUS, + SIDEBAND_DISTANCES) + rows, cols = pixels[3] + assert cols.min() == 0 + assert cols.size > 0 + + +class TestSidebandBackground: + @pytest.mark.parametrize('bayer_mask', [False, True]) + def test_constant_pedestal_is_recovered(self, bayer_mask): + search_areas = make_strip(30.0) + img, mask, _ = make_image(search_areas, pedestal=7.0, bayer_mask=bayer_mask) + pixels = get_sideband_pixels(search_areas, [search_areas[:, 0]], img.shape, SEARCH_AREA_RADIUS, + SIDEBAND_DISTANCES) + for led_pixels in pixels: + assert estimate_sideband_background(img, mask, led_pixels) == pytest.approx(7.0) + + def test_clipped_zeros_are_included(self): + # pedestal of 1 count with noise of 2 counts: clipping at 0 raises the mean in the box and in the side + # bands alike, so the side band mean (including zeros) is the matching estimator; the mean of the + # pixels > 0 is too high + from scipy.stats import norm + search_areas = make_strip(0.0) + img, mask, _ = make_image(search_areas, pedestal=1.0, noise_sigma=2.0, clip=True, bayer_mask=True) + clipped_mean = 1.0 * norm.cdf(0.5) + 2.0 * norm.pdf(0.5) # expectation of max(0, N(1, 2)) + pixels = get_sideband_pixels(search_areas, [search_areas[:, 0]], img.shape, SEARCH_AREA_RADIUS, + SIDEBAND_DISTANCES) + estimate = estimate_sideband_background(img, mask, pixels[3]) + positive_mean = np.mean(img[pixels[3]][img[pixels[3]] > 0]) + assert estimate == pytest.approx(clipped_mean, rel=0.1) + assert positive_mean > 1.3 * estimate + + def test_no_channel_pixels_gives_nan(self): + img = np.zeros((10, 10)) + assert np.isnan(estimate_sideband_background(img, None, None)) + assert np.isnan(estimate_sideband_background(img, None, (np.array([], dtype=int), np.array([], dtype=int)))) + mask = np.zeros((10, 10), dtype=bool) + assert np.isnan(estimate_sideband_background(img, mask, (np.array([1]), np.array([1])))) + + +class TestSidebandSum: + @pytest.mark.parametrize('bayer_mask', [False, True]) + @pytest.mark.parametrize('angle', [0.0, 90.0, 30.0]) + def test_sidebg_sum_is_independent_of_pedestal(self, bayer_mask, angle): + search_areas = make_strip(angle) + for pedestal in [0.0, 3.0, 20.0]: + img, mask, led_img = make_image(search_areas, pedestal=pedestal, bayer_mask=bayer_mask) + for led_id, row, col in search_areas: + led_data = led_data_for(img, mask, search_areas, led_id) + assert led_data.sidebg_sum_color_value == pytest.approx(box_sum(led_img, row, col), rel=1e-6) + num_channel_pixels = (2 * SEARCH_AREA_RADIUS) ** 2 // (2 if bayer_mask else 1) + assert led_data.sidebg_pedestal == pytest.approx(num_channel_pixels * pedestal, abs=1e-9) + + def test_optical_depth_bias_with_clipped_noise(self): + # optical depth -ln(I/I0) of a dimmed LED with pedestal, clipped noise and Bayer mask: the side band + # correction removes the bias of the plain box sum + search_areas = make_strip(0.0) + transmission = 0.01 + image_args = dict(pedestal=2.0, noise_sigma=2.0, clip=True, bayer_mask=True, amplitude=4000.0) + clean, mask, _ = make_image(search_areas, transmission=1.0, seed=1, **image_args) + smoke, _, _ = make_image(search_areas, transmission=transmission, seed=2, **image_args) + true_depth = -np.log(transmission) + depth_sum, depth_sidebg = [], [] + for led_id in search_areas[:, 0]: + ref = led_data_for(clean, mask, search_areas, led_id) + img = led_data_for(smoke, mask, search_areas, led_id) + depth_sum.append(-np.log(img.sum_color_value / ref.sum_color_value)) + depth_sidebg.append(-np.log(img.sidebg_sum_color_value / ref.sidebg_sum_color_value)) + assert abs(np.mean(depth_sidebg) - true_depth) / true_depth < 0.02 + assert (true_depth - np.mean(depth_sum)) / true_depth > 0.1 + + def test_floor_limits_the_corrected_sum(self): + assert subtract_pedestal(100.0, 99.0, 0.05) == pytest.approx(5.0) + assert subtract_pedestal(100.0, 20.0, 0.05) == pytest.approx(80.0) + assert np.isnan(subtract_pedestal(100.0, np.nan, 0.05)) + + def test_old_values_are_unchanged(self): + # the existing columns keep their values; the new ones are appended to the CSV line + search_areas = make_strip(0.0) + img, mask, _ = make_image(search_areas, pedestal=5.0, noise_sigma=1.0, bayer_mask=True) + led_data = led_data_for(img, mask, search_areas, 3) + reference = _generate_led_analysis_data(None, 1, img, False, 3, 'test_img', 0, search_areas, + SEARCH_AREA_RADIUS, fit_leds=False) + assert led_data.sum_color_value == reference.sum_color_value + assert led_data.mean_color_value == reference.mean_color_value + assert led_data.max_color_value == reference.max_color_value + assert led_data.bgsub_sum_color_value == reference.bgsub_sum_color_value + fields = led_data.get_main_data_string().split(',') + assert ','.join(fields[:6]) == ','.join(reference.get_main_data_string().split(',')[:6]) + assert float(fields[6]) == pytest.approx(led_data.sidebg_sum_color_value, rel=1e-4) + assert float(fields[7]) == pytest.approx(led_data.sidebg_pedestal, rel=1e-4) + + +class TestSidebandParameters: + def test_defaults_without_config(self): + assert get_sideband_parameters(None) == ((25, 60), 0.05, None) + + def test_defaults_for_old_config(self, tmp_path, monkeypatch): + from ledsa.core.ConfigData import ConfigData + monkeypatch.chdir(tmp_path) + (tmp_path / 'config.ini').write_text('[DEFAULT]\nimg_directory = x\n[analyse_photo]\nnum_skip_leds = 0\n', encoding='utf-8') + assert get_sideband_parameters(ConfigData()) == ((25, 60), 0.05, None) + + def test_values_from_config(self, tmp_path, monkeypatch): + from ledsa.core.ConfigData import ConfigData + monkeypatch.chdir(tmp_path) + ConfigData(load_config_file=False, sideband_distances='12 20', sideband_floor=0.1, + sideband_exclusion_margin=3) + assert get_sideband_parameters(ConfigData()) == ((12.0, 20.0), 0.1, 3) + + +class TestUnmergedLedArrays: + def test_without_merge_the_analysis_arrays_are_used(self): + conf = {'analyse_positions': {'merge_led_array_indices': 'None', 'num_arrays': '2'}} + assert _get_unmerged_led_arrays(conf, [[0, 1], [2, 3]]) == [[0, 1], [2, 3]] + + def test_with_merge_the_arrays_of_step_2_are_read(self, tmp_path, monkeypatch): + # merged arrays may contain several strips, the strip direction is estimated per array of step 2 + from ledsa.data_extraction.step_2_functions import generate_led_array_indices_files + monkeypatch.chdir(tmp_path) + generate_led_array_indices_files([[0, 1, 2], [3, 4, 5]]) + conf = {'analyse_positions': {'merge_led_array_indices': '0 1', 'num_arrays': '2'}} + led_arrays = _get_unmerged_led_arrays(conf, [[0, 1, 2, 3, 4, 5]]) + assert [sorted(led_ids.tolist()) for led_ids in led_arrays] == [[0, 1, 2], [3, 4, 5]] + + +class TestNanReferenceProperty: + @pytest.fixture + def coefficients(self, monkeypatch): + from types import SimpleNamespace + import pandas as pd + import ledsa.analysis.ExtinctionCoefficients as extinction_coefficients + from ledsa.analysis.ExtinctionCoefficientsLinear import ExtinctionCoefficientsLinear + index = pd.MultiIndex.from_product([[1, 2], [0, 1, 2]], names=['img_id', 'led_id']) + img_data = pd.DataFrame({'led_array_id': 0, 'sum_col_val': 100.0, + 'sidebg_sum_col_val': [90.0, np.nan, 90.0, 90.0, np.nan, 90.0]}, index=index) + monkeypatch.setattr(extinction_coefficients, 'read_hdf', lambda channel, path: img_data) + coefficients = ExtinctionCoefficientsLinear.__new__(ExtinctionCoefficientsLinear) + coefficients.experiment = SimpleNamespace(channel=0, path='.', led_array=0) + coefficients.average_images = False + return coefficients + + def test_nan_values_stop_the_analysis(self, coefficients): + # the solvers would fail with an unspecific error (linear) or return zeros (nonlinear) + coefficients.reference_property = 'sidebg_sum_col_val' + with pytest.raises(ValueError, match=r'2 values of sidebg_sum_col_val of led array 0 are NaN \(led ids 1\)'): + coefficients.load_img_data() + + def test_property_without_nan_is_loaded(self, coefficients): + coefficients.reference_property = 'sum_col_val' + coefficients.load_img_data() + assert coefficients.analysed_img_ids == [1, 2] + + +class TestChannelMask: + # RGGB pattern as in rawpy raw_colors_visible: 0 red, 1 and 3 green, 2 blue + FILTER_ARRAY = np.tile(np.array([[0, 1], [3, 2]]), (2, 2)) + + @pytest.mark.parametrize('channel, colors', [(0, [0]), (1, [1, 3]), (2, [2])]) + def test_bayer_channel_mask(self, channel, colors): + mask = get_bayer_channel_mask(self.FILTER_ARRAY, channel) + np.testing.assert_array_equal(mask, np.isin(self.FILTER_ARRAY, 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: (data, self.FILTER_ARRAY)) + channel_array, mask = image_reading.read_channel_data_and_mask_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)) + + 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) + 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)