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 7a47836..05389f9 100644 --- a/ledsa/analysis/ExperimentData.py +++ b/ledsa/analysis/ExperimentData.py @@ -37,6 +37,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. @@ -67,6 +76,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 @@ -113,6 +125,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 a28957f..c56eabd 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, avg_img_ids_of_ref_imgs, +from ledsa.core.file_handling import (NOISE_COLUMNS, NUM_IMGS_PER_AVERAGE, avg_img_ids_of_ref_imgs, create_analysis_infos_avg, extend_hdf, get_extinction_coefficients_file_path, read_analysis_img_ids, read_experiment_times, @@ -59,7 +59,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 @@ -76,7 +76,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 @@ -94,10 +110,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: @@ -124,7 +149,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: @@ -141,7 +170,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 @@ -200,7 +234,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}!") @@ -218,6 +259,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/file_handling.py b/ledsa/core/file_handling.py index 3ae7884..540057a 100644 --- a/ledsa/core/file_handling.py +++ b/ledsa/core/file_handling.py @@ -688,6 +688,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"] @@ -699,8 +700,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 1cd0af8..582af3c 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 @@ -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) @@ -248,25 +260,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 +300,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/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/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 1e49764..6db109d 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 __exit__(self, *args): + return False + + monkeypatch.setattr(image_reading.rawpy, 'imread', lambda filename: Opened()) + return install - 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) + +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,113 @@ 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)) + + +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_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 c0a4330..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): @@ -286,9 +289,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 +304,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)