Unclipped LED intensity, its read noise and the SNR weighting of the linear solver

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 ')

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 [])}

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:
         """