diff --git a/ledsa/core/image_reading.py b/ledsa/core/image_reading.py index b224002..fb6e1a1 100644 --- a/ledsa/core/image_reading.py +++ b/ledsa/core/image_reading.py @@ -113,6 +113,10 @@ class ImageChannels: when the object is created, and channel_data_mask_and_noise builds the arrays of a channel from it. They are built for one channel at a time, so that besides the decoded image only the channel in work is held in memory. + + The stored values of a RAW image are kept as they are, and the black level, the white level and the read + noise are applied to the pixels that are asked for. Step 3 asks only for its search areas and side bands, + a few percent of the sensor, which saves applying them to the whole image. """ def __init__(self, filename: str, rotation: str = AUTO_ROTATION): @@ -126,40 +130,82 @@ def __init__(self, filename: str, rotation: str = AUTO_ROTATION): :raises ValueError: If rotation is neither 'auto' nor one of 0, 90, 180 and 270. """ def read_img_file(img_filename, _, img_orientation=None): - # every pixel of a JPG or PNG belongs to every channel, and there is no optical black to take the - # read noise from - return _read_img_array_from_img_file(img_filename, img_orientation), None, None + return _read_img_array_from_img_file(img_filename, img_orientation) + + def read_raw_file(raw_filename, _, raw_orientation=None): + raw = _read_raw_image(raw_filename) + if raw_orientation is None: + raw_orientation = raw.orientation + # The Bayer pattern is turned with the pixels, so the filter array keeps naming the color of every pixel + return raw._replace(data=_apply_orientation(raw.data, raw_orientation), + filter_array=_apply_orientation(raw.filter_array, raw_orientation), + orientation=UPRIGHT_ORIENTATION) self.filename = filename orientation = _orientation_from_rotation(rotation) - reader = _get_reader(filename, read_img_file, _read_img_array_and_noise_from_raw_file) + reader = _get_reader(filename, read_img_file, read_raw_file) # no channel is passed, the read is the same for all of them - self._img_array, self._filter_array, self._pixel_noise = _read_or_raise(filename, reader, None, - orientation) + img = _read_or_raise(filename, reader, None, orientation) + self._raw = img if isinstance(img, RawImage) else None + self._img_array = None if self._raw is not None else img + # the calibrated pixels of the last selection, with the selection they belong to + self._calibrated = None - def channel_data_mask_and_noise(self, channel: int) -> Tuple[np.ndarray, np.ndarray, Optional[np.ndarray]]: + @property + def shape(self) -> Tuple[int, int]: + """The number of rows and columns of the upright image.""" + return (self._img_array if self._raw is None else self._raw.data).shape[:2] + + def channel_data_mask_and_noise( + self, channel: int, + pixels: Optional[Tuple[np.ndarray, np.ndarray]] = None + ) -> Tuple[np.ndarray, np.ndarray, Optional[np.ndarray]]: """ The channel array, the boolean channel mask and the read noise of every pixel of one color channel, as - read_channel_data_mask_and_noise_from_img returns them. + read_channel_data_mask_and_noise_from_img returns them, of the whole image or of selected pixels. + + The calibrated pixels of a RAW image do not depend on the channel. They are kept for the next call that + passes the same selection, the same object, as the analysis of the next channel of an image does. :param channel: The color channel to process. :type channel: int + :param pixels: Row and column indices of the selected pixels in the upright image, None for all pixels. + :type pixels: tuple[np.ndarray, np.ndarray] or None :return: The channel array, the boolean channel mask and the read noise of every pixel, all of the - same shape. The noise is None for JPG and PNG images and for sensors without an optical black - margin. + shape of the image or of the index arrays. The noise is None for JPG and PNG images and for sensors + without an optical black margin. :rtype: tuple[np.ndarray, np.ndarray, np.ndarray or None] :raises CorruptedImageError: If the image does not have the channel. """ try: - if self._filter_array is None: + if self._raw is None: + # every pixel of a JPG or PNG belongs to every channel, and there is no optical black to take the + # read noise from channel_array = self._img_array[:, :, channel] + if pixels is not None: + channel_array = channel_array[pixels] return channel_array, np.ones(channel_array.shape, dtype=bool), None - channel_mask = get_bayer_channel_mask(self._filter_array, channel) - return (np.where(channel_mask, self._img_array, 0), channel_mask, - None if self._pixel_noise is None else np.where(channel_mask, self._pixel_noise, 0.0)) + img_array, filter_array, pixel_noise = self._calibrate(pixels) + channel_mask = get_bayer_channel_mask(filter_array, channel) + return (np.where(channel_mask, img_array, 0), channel_mask, + None if pixel_noise is None else np.where(channel_mask, pixel_noise, 0.0)) except Exception as error: raise CorruptedImageError(self.filename, error) from error + def _calibrate(self, pixels: Optional[Tuple[np.ndarray, np.ndarray]] + ) -> Tuple[np.ndarray, np.ndarray, Optional[np.ndarray]]: + """ + The values with the black level subtracted and scaled to the white level, the Bayer colors and the read + noise of the selected pixels of the RAW image, as _channel_data and _pixel_read_noise yield them for the + whole image. + """ + if self._calibrated is None or self._calibrated[0] is not pixels: + raw = self._raw if pixels is None else self._raw._replace(data=self._raw.data[pixels], + filter_array=self._raw.filter_array[pixels]) + # the selection is kept with the result, so that its identity cannot pass to another selection + self._calibrated = (pixels, (_channel_data(raw), raw.filter_array, _pixel_read_noise(raw))) + return self._calibrated[1] + def get_bayer_channel_mask(filter_array: np.ndarray, channel: int) -> np.ndarray: """ @@ -394,22 +440,6 @@ def _read_img_array_from_raw_file(filename: str, channel: int, _apply_orientation(raw.filter_array, orientation)) -def _read_img_array_and_noise_from_raw_file( - filename: str, channel: int, - orientation: Optional[int] = None) -> Tuple[np.ndarray, np.ndarray, Optional[np.ndarray]]: - """ - As _read_img_array_from_raw_file, plus the read noise of every single pixel. The file is read once for - both. - """ - raw = _read_raw_image(filename) - if orientation is None: - orientation = raw.orientation - noise = _pixel_read_noise(raw) - return (_apply_orientation(_channel_data(raw), orientation), - _apply_orientation(raw.filter_array, orientation), - None if noise is None else _apply_orientation(noise, orientation)) - - def _read_img_array_from_img_file(filename: str, orientation: Optional[int] = None) -> np.ndarray: img_array = plt.imread(filename) if orientation is None: diff --git a/ledsa/data_extraction/sideband_background.py b/ledsa/data_extraction/sideband_background.py index f4220f2..4108362 100644 --- a/ledsa/data_extraction/sideband_background.py +++ b/ledsa/data_extraction/sideband_background.py @@ -169,14 +169,40 @@ def estimate_sideband_background(data: np.ndarray, channel_mask: Optional[np.nda :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: + if sideband_pixels is None: 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)) + in_channel = np.ones(values.shape, dtype=bool) if channel_mask is None else channel_mask[sideband_pixels] + return float(estimate_sideband_backgrounds(values, in_channel, [values.size])[0]) + + +def estimate_sideband_backgrounds(values: np.ndarray, channel_mask: np.ndarray, + num_pixels: Sequence[int]) -> np.ndarray: + """ + estimate_sideband_background for many LEDs at once, from the side band pixels of all of them. + + :param values: Channel values of the side band pixels of all LEDs, the pixels of one LED after those of the + previous one. + :type values: np.ndarray + :param channel_mask: Whether each of these pixels belongs to the channel. + :type channel_mask: np.ndarray + :param num_pixels: Number of side band pixels of every LED. + :type num_pixels: Sequence[int] + :return: Mean pedestal per channel pixel of every LED, NaN for an LED whose side bands contain no channel pixel. + :rtype: np.ndarray + """ + num_pixels = np.asarray(num_pixels, dtype=np.intp) + backgrounds = np.full(num_pixels.size, np.nan) + # reduceat cannot sum an empty segment, and the segments left are contiguous without the empty ones + has_pixels = num_pixels > 0 + if not has_pixels.any(): + return backgrounds + starts = (np.cumsum(num_pixels) - num_pixels)[has_pixels] + sums = np.add.reduceat(np.where(channel_mask, values, 0), starts, dtype=np.float64) + counts = np.add.reduceat(channel_mask, starts, dtype=np.intp) + with np.errstate(invalid='ignore'): + backgrounds[has_pixels] = sums / counts # 0 / 0 is NaN + return backgrounds def subtract_pedestal(sum_color_value: float, pedestal: float, floor: float = DEFAULT_BACKGROUND_FLOOR) -> float: diff --git a/ledsa/data_extraction/step_3_functions.py b/ledsa/data_extraction/step_3_functions.py index 74f0237..4bf9138 100644 --- a/ledsa/data_extraction/step_3_functions.py +++ b/ledsa/data_extraction/step_3_functions.py @@ -2,7 +2,7 @@ import re import time from importlib.metadata import version -from typing import List, Optional, Tuple +from typing import List, NamedTuple, Optional, Sequence, Tuple import numpy as np import scipy.optimize @@ -14,8 +14,8 @@ 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, - estimate_sideband_background, get_sideband_pixels, - subtract_pedestal) + estimate_sideband_background, estimate_sideband_backgrounds, + get_sideband_pixels, subtract_pedestal) def generate_analysis_data(img_filename: str, channel: int, search_areas: np.ndarray, line_indices: List[List[int]], @@ -49,33 +49,139 @@ def generate_analysis_data(img_filename: str, channel: int, search_areas: np.nda if img_channels is None: file_path = os.path.join(conf['DEFAULT']['img_directory'], img_filename) img_channels = ImageChannels(file_path, rotation=conf.get('DEFAULT', 'img_rotation')) - data, channel_mask, pixel_noise = img_channels.channel_data_mask_and_noise(channel) search_area_radius = int(conf['find_search_areas']['search_area_radius']) - img_analysis_data = [] if debug: + data, _, _ = img_channels.channel_data_mask_and_noise(channel) analysis_res = _generate_led_analysis_data(conf, channel, data, debug, debug_led, img_filename, 0, search_areas, search_area_radius, fit_leds) return analysis_res sideband_distances, background_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, - channel_mask=channel_mask, - pixel_noise=pixel_noise, - sideband_pixels=sideband_pixels[iled], - background_floor=background_floor) - img_analysis_data.append(led_analysis_data) + sideband_pixels = get_sideband_pixels(search_areas, _get_unmerged_led_arrays(conf, line_indices), + img_channels.shape, search_area_radius, sideband_distances, + sideband_exclusion_margin) + + num_skip_leds = int(conf['analyse_photo']['num_skip_leds']) + leds = [(int(iled), led_array_idx) for led_array_idx, led_array in enumerate(line_indices) + for iled in led_array if iled % (num_skip_leds + 1) == 0] + led_ids = np.array([iled for iled, _ in leds], dtype=int) + led_array_ids = np.array([led_array_idx for _, led_array_idx in leds], dtype=int) + + # only the pixels of the search areas and side bands are read, which are a few percent of the image + analysis_pixels = _get_analysis_pixels(search_areas, led_ids, sideband_pixels, img_channels.shape, + search_area_radius) + data, channel_mask, pixel_noise = img_channels.channel_data_mask_and_noise(channel, analysis_pixels.pixels) + sidebands = np.s_[analysis_pixels.num_box_pixels:] + sideband_backgrounds = estimate_sideband_backgrounds(data[sidebands], channel_mask[sidebands], + analysis_pixels.num_sideband_pixels) + + img_analysis_data = [None] * len(leds) + for positions, shape, start in analysis_pixels.box_groups: + boxes = np.s_[start:start + positions.size * shape[0] * shape[1]] + stack_shape = (positions.size,) + shape + group_data = _analyse_search_areas( + conf, channel, img_filename, led_ids[positions], led_array_ids[positions], + analysis_pixels.centers[positions], search_area_radius, fit_leds, + data[boxes].reshape(stack_shape), channel_mask[boxes].reshape(stack_shape), + None if pixel_noise is None else pixel_noise[boxes].reshape(stack_shape), + sideband_backgrounds[positions], background_floor) + for position, led_analysis_data in zip(positions, group_data): + img_analysis_data[position] = led_analysis_data return img_analysis_data +class AnalysisPixels(NamedTuple): + """ + The pixels step 3 reads from an image, in one selection so that the image calibrates them once for all + channels: the search areas of the analysed LEDs, stacked by their shape, followed by their side bands. + """ + + pixels: Tuple[np.ndarray, np.ndarray] # row and column indices + centers: np.ndarray # center (row, column) of the search area of every analysed LED + # per shape of the search areas: the positions of their LEDs in the analysis, the shape and the index of the + # first pixel of the stack in pixels + box_groups: List[Tuple[np.ndarray, Tuple[int, int], int]] + num_box_pixels: int # the side band pixels follow the search areas, LED after LED + num_sideband_pixels: np.ndarray # number of side band pixels of every analysed LED + + +# pixels of the analysis, cached per process because they only depend on the LED positions +_analysis_pixels_cache = {} + + +def _get_analysis_pixels(search_areas: np.ndarray, led_ids: np.ndarray, + sideband_pixels: List[Optional[Tuple[np.ndarray, np.ndarray]]], + img_shape: Tuple[int, int], search_area_radius: int) -> AnalysisPixels: + """ + Cached _build_analysis_pixels. The same object is returned for every image and channel, which lets the image + keep the pixels it calibrated for the next channel. + + :return: See _build_analysis_pixels. + """ + key = (search_areas.tobytes(), search_areas.shape, led_ids.tobytes(), tuple(img_shape), search_area_radius) + # the side band pixels come from a cache of their own, which builds a new list when its key changes; the + # list is kept with the entry, so that its identity cannot pass to another one + if key not in _analysis_pixels_cache or _analysis_pixels_cache[key][0] is not sideband_pixels: + _analysis_pixels_cache.clear() + _analysis_pixels_cache[key] = (sideband_pixels, _build_analysis_pixels( + search_areas, led_ids, sideband_pixels, img_shape, search_area_radius)) + return _analysis_pixels_cache[key][1] + + +def _build_analysis_pixels(search_areas: np.ndarray, led_ids: np.ndarray, + sideband_pixels: List[Optional[Tuple[np.ndarray, np.ndarray]]], + img_shape: Tuple[int, int], search_area_radius: int) -> AnalysisPixels: + """ + Collect the pixels of the search areas and side bands of the analysed LEDs. + + The search area is cut from the image as data[row - R:row + R, col - R:col + R], so an area at the border + of the image is cut short as the slice is. Its shape differs from that of the others, and the search areas + are therefore stacked by their shape. + + :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_ids: Ids of the analysed LEDs, in the order of the analysis. + :type led_ids: np.ndarray + :param sideband_pixels: Per LED id the row and column indices of its side band pixels, see get_sideband_pixels. + :type sideband_pixels: list[tuple[np.ndarray, np.ndarray] or None] + :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 + :return: The pixels and where the analysis finds the search area and the side bands of every LED in them. + :rtype: AnalysisPixels + """ + centers = search_areas[led_ids, 1:3].astype(int) + row_ranges = [range(*slice(row - search_area_radius, row + search_area_radius).indices(img_shape[0])) + for row in centers[:, 0]] + col_ranges = [range(*slice(col - search_area_radius, col + search_area_radius).indices(img_shape[1])) + for col in centers[:, 1]] + positions_by_shape = {} + for position, (row_range, col_range) in enumerate(zip(row_ranges, col_ranges)): + positions_by_shape.setdefault((len(row_range), len(col_range)), []).append(position) + + rows, cols, box_groups, start = [np.zeros(0, dtype=np.intp)], [np.zeros(0, dtype=np.intp)], [], 0 + for shape, positions in positions_by_shape.items(): + positions = np.array(positions, dtype=int) + stack_shape = (positions.size,) + shape + first_rows = np.array([row_ranges[position].start for position in positions], dtype=np.intp) + first_cols = np.array([col_ranges[position].start for position in positions], dtype=np.intp) + rows.append(np.broadcast_to(first_rows[:, None, None] + np.arange(shape[0])[:, None], stack_shape).ravel()) + cols.append(np.broadcast_to(first_cols[:, None, None] + np.arange(shape[1]), stack_shape).ravel()) + box_groups.append((positions, shape, start)) + start += positions.size * shape[0] * shape[1] + + num_sideband_pixels = np.zeros(led_ids.size, dtype=np.intp) + for position, led_id in enumerate(led_ids): + if sideband_pixels[led_id] is not None: + rows.append(sideband_pixels[led_id][0]) + cols.append(sideband_pixels[led_id][1]) + num_sideband_pixels[position] = sideband_pixels[led_id][0].size + return AnalysisPixels((np.concatenate(rows), np.concatenate(cols)), centers, box_groups, start, + num_sideband_pixels) + + 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 @@ -185,7 +291,8 @@ def _generate_led_analysis_data(conf: ConfigData, channel: int, data: np.ndarray sideband_pixels: Optional[Tuple[np.ndarray, np.ndarray]] = None, background_floor: float = DEFAULT_BACKGROUND_FLOOR) -> LEDAnalysisData: """ - Generate analysis data for a specific LED. + Generate analysis data for a specific LED from the arrays of the whole image, as _analyse_search_areas does + for many. :param conf: Configuration data. :type conf: ConfigData @@ -193,7 +300,7 @@ def _generate_led_analysis_data(conf: ConfigData, channel: int, data: np.ndarray :type channel: int :param data: Array representing the image data. :type data: np.ndarray - :param debug: If True, enters debug mode. + :param debug: If True, the fit results are returned instead of the analysis data. :type debug: bool :param iled: ID of the LED for which analysis should be generated. :type iled: int @@ -220,124 +327,155 @@ def _generate_led_analysis_data(conf: ConfigData, channel: int, data: np.ndarray :return: Analysis data for the LED. :rtype: LEDAnalysisData """ - led_data = LEDAnalysisData(iled, led_array_idx, fit_leds) - center_search_area_x = int(search_areas[iled, 1]) - center_search_area_y = int(search_areas[iled, 2]) - search_area = np.index_exp[center_search_area_x - search_area_radius: - center_search_area_x + search_area_radius, - center_search_area_y - search_area_radius: - center_search_area_y + search_area_radius] + center = np.array([[int(search_areas[iled, 1]), int(search_areas[iled, 2])]]) + search_area = np.index_exp[center[0, 0] - search_area_radius:center[0, 0] + search_area_radius, + center[0, 1] - search_area_radius:center[0, 1] + search_area_radius] + box_mask = np.ones(data[search_area].shape, dtype=bool) if channel_mask is None else channel_mask[search_area] + led_data, = _analyse_search_areas(conf, channel, img_filename, [iled], [led_array_idx], center, + search_area_radius, fit_leds, data[search_area][None], box_mask[None], + None if pixel_noise is None else pixel_noise[search_area][None], + [estimate_sideband_background(data, channel_mask, sideband_pixels)], + background_floor) + if debug and fit_leds: + return led_data.fit_results.x + return led_data - if fit_leds: - start_time = time.process_time() - led_data.fit_results, mesh = _fit_model_to_led(data[search_area]) - end_time = time.process_time() - led_data.fit_time = end_time - start_time - led_data.led_center_x = led_data.fit_results.x[0] + center_search_area_x - search_area_radius - led_data.led_center_y = led_data.fit_results.x[1] + center_search_area_y - search_area_radius - if debug: - return led_data.fit_results.x - if not led_data.fit_results.success: # A > 255 or A < 0: - _log_warnings(img_filename, channel, led_data, center_search_area_x, center_search_area_y, - data[search_area].shape, search_area_radius, conf) - - led_data.mean_color_value = np.mean(data[search_area]) - led_data.sum_color_value = np.sum(data[search_area]) - led_data.max_color_value = np.amax(data[search_area]) - box_mask = None if channel_mask is None else channel_mask[search_area] + +def _analyse_search_areas(conf: ConfigData, channel: int, img_filename: str, led_ids: Sequence[int], + led_array_ids: Sequence[int], centers: np.ndarray, search_area_radius: int, fit_leds: bool, + boxes: np.ndarray, box_masks: np.ndarray, box_noise: Optional[np.ndarray], + sideband_backgrounds: Sequence[float], + background_floor: float = DEFAULT_BACKGROUND_FLOOR) -> List[LEDAnalysisData]: + """ + Generate the analysis data of the LEDs of a stack of search areas of one shape. + + Every quantity is taken along the axes of the stack, so that an image of a thousand LEDs needs a few numpy + calls per quantity instead of a few per LED and quantity. Only the fit is made LED by LED. + + :param conf: Configuration data, used to log irregular fits. + :type conf: ConfigData + :param channel: Color channel being analysed. + :type channel: int + :param img_filename: Name of the image. + :type img_filename: str + :param led_ids: ID of the LED of every search area. + :type led_ids: Sequence[int] + :param led_array_ids: Index of the LED array of every LED. + :type led_array_ids: Sequence[int] + :param centers: Center (row, column) of every search area in the image. + :type centers: np.ndarray + :param search_area_radius: Radius of the search areas. + :type search_area_radius: int + :param fit_leds: If True, every LED is fitted to a model function. + :type fit_leds: bool + :param boxes: Channel values of the search areas, 0 outside the channel, of shape (number of LEDs, rows, + columns). + :type boxes: np.ndarray + :param box_masks: Whether each pixel of the search areas belongs to the channel, of the shape of boxes. + :type box_masks: np.ndarray + :param box_noise: Read noise of each pixel of the search areas, 0 outside the channel, None if the image has + no optical black margin. If None, the noise of the sums stays NaN. + :type box_noise: np.ndarray or None + :param sideband_backgrounds: Mean pedestal per channel pixel in the side bands of every LED, see + estimate_sideband_background. + :type sideband_backgrounds: Sequence[float] + :param background_floor: Lower limit of the side band and ring corrected sums as fraction of the box sum. + :type background_floor: float + :return: Analysis data for every LED. + :rtype: List[LEDAnalysisData] + """ + means = np.mean(boxes, axis=(1, 2)) + sums = np.sum(boxes, axis=(1, 2)) + maxima = np.amax(boxes, axis=(1, 2)) # The median of the border ring, subtracted from every channel pixel of the box. It used to be # clipped at zero pixel by pixel, which puts back the offset that the unclipped data remove. - background = _estimate_local_background(data[search_area], box_mask) - pedestal_per_pixel = background if box_mask is None else background * box_mask - led_data.bgsub_sum_color_value = float(np.sum(data[search_area].astype(np.float64) - - pedestal_per_pixel)) + ring_medians, ring_means = _estimate_ring_backgrounds(boxes, box_masks) + bgsub_sums = np.sum(boxes.astype(np.float64) - ring_medians[:, None, None] * box_masks, axis=(1, 2)) - if pixel_noise is not None: - # the noise is 0 outside the channel, so the variance of the sum needs no mask - led_data.noise_sum_color_value = float( - np.sqrt(np.sum(pixel_noise[search_area].astype(np.float64) ** 2))) + # the noise is 0 outside the channel, so the variance of the sum needs no mask + noise_sums = None if box_noise is None else np.sqrt(np.sum(box_noise.astype(np.float64) ** 2, axis=(1, 2))) # the pedestals below are subtracted as N * b with N the number of channel pixels in the box - num_channel_pixels = data[search_area].size if box_mask is None else int(np.count_nonzero(box_mask)) + num_channel_pixels = np.count_nonzero(box_masks, axis=(1, 2)) - # 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, background_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, - background_floor) - - return led_data + img_analysis_data = [] + for i, led_id in enumerate(led_ids): + led_data = LEDAnalysisData(led_id, int(led_array_ids[i]), fit_leds) + center_search_area_x, center_search_area_y = int(centers[i, 0]), int(centers[i, 1]) + if fit_leds: + start_time = time.process_time() + led_data.fit_results, mesh = _fit_model_to_led(boxes[i]) + end_time = time.process_time() + led_data.fit_time = end_time - start_time + led_data.led_center_x = led_data.fit_results.x[0] + center_search_area_x - search_area_radius + led_data.led_center_y = led_data.fit_results.x[1] + center_search_area_y - search_area_radius + if not led_data.fit_results.success: # A > 255 or A < 0: + _log_warnings(img_filename, channel, led_data, center_search_area_x, center_search_area_y, + boxes[i].shape, search_area_radius, conf) + + led_data.mean_color_value = means[i] + led_data.sum_color_value = sums[i] + led_data.max_color_value = maxima[i] + led_data.bgsub_sum_color_value = float(bgsub_sums[i]) + if noise_sums is not None: + led_data.noise_sum_color_value = float(noise_sums[i]) + + # 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 + led_data.bgsub2_sum_color_value = subtract_pedestal(float(sums[i]), num_channel_pixels[i] * ring_means[i], + background_floor) + + # b = mean over the channel pixels in the side bands next to the LED strip + led_data.sidebg_pedestal = float(num_channel_pixels[i] * sideband_backgrounds[i]) + led_data.sidebg_sum_color_value = subtract_pedestal(float(sums[i]), led_data.sidebg_pedestal, + background_floor) + img_analysis_data.append(led_data) + return img_analysis_data -def _estimate_local_background(search_area_data: np.ndarray, channel_mask: Optional[np.ndarray] = None, - border_width: int = 2) -> float: +def _estimate_ring_backgrounds(boxes: np.ndarray, box_masks: Optional[np.ndarray] = None, + border_width: int = 2) -> Tuple[np.ndarray, np.ndarray]: """ - Estimate the local background of a search area from the median of its border pixels. + Estimate the local background of every search area of a stack from the channel pixels of its border ring, + as their median and as their mean. The LED sits at the center of the search area, so the outer border ring samples the scene background (stray light, smoke path radiance). Pixels of foreign Bayer colors are excluded through the channel mask and no longer by their value: the data are not clipped at zero any more, so a pixel of the channel may well be 0 or negative, and dropping those would raise the median. - :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: Estimated background value per pixel. 0 if no valid border pixels exist. - :rtype: float - """ - if search_area_data.shape[0] <= 2 * border_width or search_area_data.shape[1] <= 2 * border_width: - return 0.0 - 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 - border_pixels = search_area_data[border_mask] - if border_pixels.size == 0: - return 0.0 - 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 + The median is subtracted from the box pixel by pixel (bgsub), the mean as N * b from the box sum (N channel + pixels in the box, bgsub2), which is linear in the pixel values like the side band correction. Note that the + border ring lies in the halo of the LED for small search areas. + + :param boxes: Search areas of one shape, stacked along the first axis. + :type boxes: np.ndarray + :param box_masks: Boolean mask of the channel pixels of the search areas, None if all pixels belong to the channel. - :type channel_mask: np.ndarray or None + :type box_masks: 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)) + :return: The median of every search area, 0 if its ring contains no channel pixel, and the mean, NaN then. + :rtype: tuple[np.ndarray, np.ndarray] + """ + medians = np.zeros(boxes.shape[0]) + means = np.full(boxes.shape[0], np.nan) + if boxes.shape[1] <= 2 * border_width or boxes.shape[2] <= 2 * border_width: + return medians, means + in_ring = np.ones(boxes.shape, dtype=bool) + in_ring[:, border_width:-border_width, border_width:-border_width] = False + if box_masks is not None: + in_ring &= box_masks + num_ring_pixels = np.count_nonzero(in_ring, axis=(1, 2)) + # The rings of a Bayer pattern hold the same number of channel pixels wherever the search area lies, so + # the search areas are grouped by that number and the estimates are taken along the axis of a group. + for num in np.unique(num_ring_pixels[num_ring_pixels > 0]): + group = num_ring_pixels == num + ring_pixels = boxes[group][in_ring[group]].reshape(-1, num) + medians[group] = np.median(ring_pixels, axis=1) + means[group] = np.mean(ring_pixels, axis=1, dtype=np.float64) + return medians, means def _save_results_in_file(channel: int, img_data: LEDAnalysisData, img_filename: str, img_id: str, img_infos: np.ndarray, basename: str) -> None: diff --git a/ledsa/tests/UnitTests/test_background_subtraction.py b/ledsa/tests/UnitTests/test_background_subtraction.py index 7c19cb7..3c7b461 100644 --- a/ledsa/tests/UnitTests/test_background_subtraction.py +++ b/ledsa/tests/UnitTests/test_background_subtraction.py @@ -1,7 +1,7 @@ import numpy as np import pytest -from ledsa.data_extraction.step_3_functions import _estimate_local_background, _generate_led_analysis_data +from ledsa.data_extraction.step_3_functions import _estimate_ring_backgrounds, _generate_led_analysis_data SEARCH_AREA_RADIUS = 10 LED_CENTER = (30, 30) @@ -9,6 +9,12 @@ LED_SIGMA = 2.0 +def ring_median(search_area, channel_mask=None): + """The median of the border ring of a single search area.""" + medians, _ = _estimate_ring_backgrounds(search_area[None], None if channel_mask is None else channel_mask[None]) + return medians[0] + + def make_synthetic_image(background=0.0, bayer_mask=False): """Create a 60x60 image with a Gaussian LED spot at LED_CENTER plus a constant background. @@ -43,10 +49,10 @@ def extract_led_data(img, channel_mask=None): SEARCH_AREA_RADIUS, fit_leds=False, channel_mask=channel_mask) -class TestEstimateLocalBackground: +class TestRingMedian: def test_constant_image_returns_background(self): img = np.full((20, 20), 7.5) - assert _estimate_local_background(img) == pytest.approx(7.5) + assert ring_median(img) == pytest.approx(7.5) def test_foreign_pixels_are_excluded_by_the_mask_and_not_by_their_value(self): # the data are no longer clipped at zero, so a channel pixel may be 0 or negative and only the @@ -55,21 +61,21 @@ def test_foreign_pixels_are_excluded_by_the_mask_and_not_by_their_value(self): channel_mask = np.zeros(img.shape, dtype=bool) channel_mask[::2, ::2] = True img[~channel_mask] = 0.0 - assert _estimate_local_background(img, channel_mask) == pytest.approx(7.5) - assert _estimate_local_background(img) == pytest.approx(0.0) # without the mask they all count + assert ring_median(img, channel_mask) == pytest.approx(7.5) + assert ring_median(img) == pytest.approx(0.0) # without the mask they all count def test_channel_pixels_at_or_below_zero_are_kept(self): img = np.full((20, 20), -1.5) - assert _estimate_local_background(img) == pytest.approx(-1.5) + assert ring_median(img) == pytest.approx(-1.5) def test_all_zero_border_returns_zero(self): img = np.zeros((20, 20)) img[10, 10] = 100.0 - assert _estimate_local_background(img) == 0.0 + assert ring_median(img) == 0.0 def test_too_small_search_area_returns_zero(self): img = np.full((4, 4), 7.5) - assert _estimate_local_background(img) == 0.0 + assert ring_median(img) == 0.0 class TestBackgroundSubtractedSum: diff --git a/ledsa/tests/UnitTests/test_black_level.py b/ledsa/tests/UnitTests/test_black_level.py index ad13e7f..da2e37c 100644 --- a/ledsa/tests/UnitTests/test_black_level.py +++ b/ledsa/tests/UnitTests/test_black_level.py @@ -4,9 +4,8 @@ import pytest import ledsa.core.image_reading as image_reading -from ledsa.core.image_reading import (MIN_OPTICAL_BLACK_WIDTH, UPRIGHT_ORIENTATION, _channel_data, - _optical_black_statistics, _pixel_read_noise, - _read_img_array_and_noise_from_raw_file, _read_raw_image) +from ledsa.core.image_reading import (MIN_OPTICAL_BLACK_WIDTH, ImageChannels, _channel_data, + _optical_black_statistics, _pixel_read_noise, _read_raw_image) # RGGB pattern as in rawpy raw_colors: 0 red, 1 and 3 green, 2 blue BAYER_TILE = np.array([[0, 1], [3, 2]]) @@ -141,16 +140,16 @@ def test_there_is_no_noise_without_an_optical_black_margin(self, fake_imread): class TestTheReaderStepThreeUses: def test_the_orientation_reaches_the_data_the_filter_and_the_noise(self, fake_imread): fake_imread(make_raw(MIN_OPTICAL_BLACK_WIDTH, noise=4.0, visible=5000.0)) - upright = _read_img_array_and_noise_from_raw_file('img.CR3', 0, UPRIGHT_ORIENTATION) - turned = _read_img_array_and_noise_from_raw_file('img.CR3', 0, 6) + upright = ImageChannels('img.CR3', rotation=0).channel_data_mask_and_noise(0) + turned = ImageChannels('img.CR3', rotation=90).channel_data_mask_and_noise(0) for upright_array, turned_array in zip(upright, turned): np.testing.assert_array_equal(turned_array, np.rot90(upright_array, -1)) def test_without_a_margin_only_the_noise_is_missing(self, fake_imread): fake_imread(make_raw(MIN_OPTICAL_BLACK_WIDTH - 1, visible=5000.0)) - data, filter_array, noise = _read_img_array_and_noise_from_raw_file('img.CR3', 0) + data, channel_mask, noise = ImageChannels('img.CR3').channel_data_mask_and_noise(0) assert noise is None - assert data.shape == filter_array.shape and np.isfinite(data).all() + assert data.shape == channel_mask.shape and np.isfinite(data).all() class TestSumsOfAnLed: diff --git a/ledsa/tests/UnitTests/test_image_channels.py b/ledsa/tests/UnitTests/test_image_channels.py index 29de7ae..9ad47ea 100644 --- a/ledsa/tests/UnitTests/test_image_channels.py +++ b/ledsa/tests/UnitTests/test_image_channels.py @@ -103,6 +103,51 @@ def test_unsupported_format_raises(self): ImageChannels('img.xyz') +class TestSelectedPixels: + @pytest.mark.parametrize('orientation', [image_reading.UPRIGHT_ORIENTATION, SIDEWAYS]) + def test_raw_pixels_are_those_of_the_whole_image(self, raw_reads, orientation): + raw_reads(orientation) + img_channels = ImageChannels('img.CR3') + pixels = (np.array([0, 3, 1, 1]), np.array([2, 0, 3, 3])) # in the upright image, one of them twice + for channel in CHANNEL_COLORS: + for selected, whole in zip(img_channels.channel_data_mask_and_noise(channel, pixels), + img_channels.channel_data_mask_and_noise(channel)): + np.testing.assert_array_equal(selected, whole[pixels]) + + def test_png_pixels_are_those_of_the_color_planes(self, tmp_path): + path = str(tmp_path / 'img.png') + plt.imsave(path, np.random.default_rng(0).random((6, 8, 3))) + img_channels = ImageChannels(path) + pixels = (np.array([[0, 5], [2, 2]]), np.array([[7, 0], [3, 4]])) # index arrays of any shape + for channel in CHANNEL_COLORS: + channel_array, channel_mask, pixel_noise = img_channels.channel_data_mask_and_noise(channel, pixels) + np.testing.assert_array_equal(channel_array, img_channels.channel_data_mask_and_noise(channel)[0][pixels]) + assert channel_mask.shape == channel_array.shape == (2, 2) and channel_mask.all() + assert pixel_noise is None + + def test_the_calibration_is_kept_for_the_same_selection_only(self, raw_reads, monkeypatch): + raw_reads() + calibrated = [] + channel_data = image_reading._channel_data + monkeypatch.setattr(image_reading, '_channel_data', + lambda raw: calibrated.append(raw.data.shape) or channel_data(raw)) + img_channels = ImageChannels('img.CR3') + pixels = (np.array([0, 1]), np.array([2, 3])) + for channel in CHANNEL_COLORS: + img_channels.channel_data_mask_and_noise(channel, pixels) + assert calibrated == [(2,)] + # another selection, even of the same pixels, and the whole image are calibrated anew + img_channels.channel_data_mask_and_noise(0, (pixels[0].copy(), pixels[1].copy())) + img_channels.channel_data_mask_and_noise(0) + img_channels.channel_data_mask_and_noise(1) + assert calibrated == [(2,), (2,), FILTER_ARRAY.shape] + + def test_shape_is_that_of_the_upright_image(self, raw_reads): + raw_reads(SIDEWAYS) + assert ImageChannels('img.CR3').shape == FILTER_ARRAY.shape[::-1] + assert ImageChannels('img.CR3', rotation=0).shape == FILTER_ARRAY.shape + + class TestStep3ReadsEveryImageOnce: @pytest.fixture def extractor(self, tmp_path, monkeypatch): diff --git a/ledsa/tests/UnitTests/test_ring_background_mean.py b/ledsa/tests/UnitTests/test_ring_background_mean.py index c467560..c2de990 100644 --- a/ledsa/tests/UnitTests/test_ring_background_mean.py +++ b/ledsa/tests/UnitTests/test_ring_background_mean.py @@ -1,13 +1,24 @@ 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) +from ledsa.data_extraction.step_3_functions import _estimate_ring_backgrounds, _generate_led_analysis_data SEARCH_AREA_RADIUS = 8 CENTER = (40, 40) +def ring_median(search_area, channel_mask=None): + """The median of the border ring of a single search area.""" + medians, _ = _estimate_ring_backgrounds(search_area[None], None if channel_mask is None else channel_mask[None]) + return medians[0] + + +def ring_mean(search_area, channel_mask=None): + """The mean of the border ring of a single search area.""" + _, means = _estimate_ring_backgrounds(search_area[None], None if channel_mask is None else channel_mask[None]) + return means[0] + + 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.""" @@ -34,24 +45,41 @@ def box(array): class TestRingBackgroundMean: def test_constant_image(self): - assert _estimate_local_background_mean(np.full((16, 16), 7.5)) == pytest.approx(7.5) + assert ring_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) + assert ring_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 ring_mean(data) == pytest.approx(data[ring].mean()) # both estimators keep the zeros now; only the Bayer mask decides which pixels count - assert _estimate_local_background(data) == pytest.approx(data[ring].mean()) + assert ring_median(data) == pytest.approx(data[ring].mean()) 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))) + assert np.isnan(ring_mean(np.ones((4, 4)))) + assert np.isnan(ring_mean(np.ones((16, 16)), np.zeros((16, 16), dtype=bool))) + + +class TestRingBackgroundsOfAStack: + def test_search_areas_with_different_numbers_of_ring_pixels(self): + rng = np.random.default_rng(2) + boxes = rng.normal(5.0, 2.0, (6, 12, 12)).astype(np.float32) + masks = rng.random(boxes.shape) < 0.5 + masks[4] = False # no channel pixel at all + ring = np.ones((12, 12), dtype=bool) + ring[2:-2, 2:-2] = False + medians, means = _estimate_ring_backgrounds(boxes, masks) + for search_area, mask, median, mean in zip(boxes, masks, medians, means): + ring_pixels = search_area[ring & mask] + assert median == (np.median(ring_pixels) if ring_pixels.size else 0.0) + assert mean == pytest.approx(np.mean(ring_pixels, dtype=np.float64) if ring_pixels.size else np.nan, + nan_ok=True) + assert medians[4] == 0.0 and np.isnan(means[4]) class TestBgsub2Sum: @@ -60,7 +88,7 @@ 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)) + expected = box(img).sum() - num_channel_pixels * ring_mean(box(img), box(mask)) assert led_data.bgsub2_sum_color_value == pytest.approx(expected) def _box_sum_below_zero(self): @@ -69,8 +97,7 @@ def _box_sum_below_zero(self): 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 num_channel_pixels = np.count_nonzero(box(mask)) - corrected = box(img).sum() - num_channel_pixels * _estimate_local_background_mean(box(img), - box(mask)) + corrected = box(img).sum() - num_channel_pixels * ring_mean(box(img), box(mask)) assert corrected < 0 return img, mask, corrected diff --git a/ledsa/tests/UnitTests/test_sideband_background.py b/ledsa/tests/UnitTests/test_sideband_background.py index 8a0030e..35e3597 100644 --- a/ledsa/tests/UnitTests/test_sideband_background.py +++ b/ledsa/tests/UnitTests/test_sideband_background.py @@ -2,8 +2,8 @@ import pytest from ledsa.data_extraction.sideband_background import (build_sideband_pixels, estimate_sideband_background, - estimate_strip_directions, get_sideband_pixels, - subtract_pedestal) + estimate_sideband_backgrounds, 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, @@ -161,6 +161,38 @@ def test_no_channel_pixels_gives_nan(self): assert np.isnan(estimate_sideband_background(img, mask, (np.array([1]), np.array([1])))) +class TestSidebandBackgroundsOfManyLeds: + def test_every_led_gets_the_mean_of_its_own_channel_pixels(self): + # the second and the last LED have no side band pixels, the third none of the channel + num_pixels = [5, 0, 3, 4, 0] + rng = np.random.default_rng(1) + values = rng.normal(10.0, 3.0, sum(num_pixels)).astype(np.float32) + mask = rng.random(values.size) < 0.6 + mask[[0, 8]] = True + mask[5:8] = False + backgrounds = estimate_sideband_backgrounds(values, mask, num_pixels) + starts = np.cumsum(num_pixels) - num_pixels + for background, start, num in zip(backgrounds, starts, num_pixels): + led_values, led_mask = values[start:start + num], mask[start:start + num] + expected = np.mean(led_values[led_mask], dtype=np.float64) if led_mask.any() else np.nan + assert background == pytest.approx(expected, nan_ok=True) + assert np.isnan(backgrounds[[1, 2, 4]]).all() + + def test_no_side_band_pixels_at_all(self): + assert np.isnan(estimate_sideband_backgrounds(np.zeros(0), np.zeros(0, dtype=bool), [0, 0])).all() + + def test_one_led_is_estimate_sideband_background(self): + search_areas = make_strip(30.0) + img, mask, _ = make_image(search_areas, pedestal=3.0, noise_sigma=1.0, bayer_mask=True) + pixels = get_sideband_pixels(search_areas, [search_areas[:, 0]], img.shape, SEARCH_AREA_RADIUS, + SIDEBAND_DISTANCES) + backgrounds = estimate_sideband_backgrounds(np.concatenate([img[led_pixels] for led_pixels in pixels]), + np.concatenate([mask[led_pixels] for led_pixels in pixels]), + [led_pixels[0].size for led_pixels in pixels]) + assert backgrounds == pytest.approx([estimate_sideband_background(img, mask, led_pixels) + for led_pixels in pixels]) + + class TestSidebandSum: @pytest.mark.parametrize('bayer_mask', [False, True]) @pytest.mark.parametrize('angle', [0.0, 90.0, 30.0]) diff --git a/ledsa/tests/UnitTests/test_stacked_search_areas.py b/ledsa/tests/UnitTests/test_stacked_search_areas.py new file mode 100644 index 0000000..35b41b6 --- /dev/null +++ b/ledsa/tests/UnitTests/test_stacked_search_areas.py @@ -0,0 +1,126 @@ +"""Unit tests for the analysis of all search areas of an image at once. + +Step 3 reads only the pixels of the search areas and side bands from an image, calibrates them once for all +channels and takes every quantity along the axes of a stack of search areas. The results have to be those of +the definitions applied LED by LED to the arrays of the whole image. +""" +import configparser + +import numpy as np +import pytest + +from ledsa.core import image_reading +from ledsa.core.image_reading import ImageChannels, RawImage +from ledsa.data_extraction.sideband_background import get_sideband_pixels +from ledsa.data_extraction.step_3_functions import (_get_analysis_pixels, generate_analysis_data, + get_sideband_parameters) + +SEARCH_AREA_RADIUS = 4 +IMG_SHAPE = (40, 48) +# 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]]), (IMG_SHAPE[0] // 2, IMG_SHAPE[1] // 2)) +SIDEWAYS = 6 +# A vertical strip and an LED alone in its array, whose strip direction and therefore side bands are unknown. +# The search area of LED 4 is cut short by the lower border of the upright image, that of LED 5 by the right +# border of the image turned sideways. +SEARCH_AREAS = np.array([[0, 8, 20], [1, 16, 20], [2, 24, 20], [3, 32, 20], [4, 38, 20], [5, 20, 38]]) +LINE_INDICES = [np.array([0, 1, 2, 3, 4]), np.array([5])] +QUANTITIES = ['sum_color_value', 'mean_color_value', 'max_color_value', 'bgsub_sum_color_value', + 'sidebg_sum_color_value', 'sidebg_pedestal', 'bgsub2_sum_color_value', 'noise_sum_color_value'] + + +def make_config(): + config = configparser.ConfigParser() + config['DEFAULT'] = {'img_directory': '.', 'img_rotation': 'auto'} + config['find_search_areas'] = {'search_area_radius': str(SEARCH_AREA_RADIUS)} + config['analyse_photo'] = {'num_skip_leds': '0', 'sideband_distances': '6 12', 'background_floor': 'None', + 'sideband_exclusion_margin': 'None'} + config['analyse_positions'] = {'merge_led_array_indices': 'None'} + return config + + +def sideband_pixels_of(img_channels, config): + distances, _, margin = get_sideband_parameters(config) + return get_sideband_pixels(SEARCH_AREAS, LINE_INDICES, img_channels.shape, SEARCH_AREA_RADIUS, distances, margin) + + +@pytest.fixture(params=[image_reading.UPRIGHT_ORIENTATION, SIDEWAYS]) +def img_channels(request, monkeypatch): + """A RAW image with a black level and a read noise of its own for every Bayer color.""" + stored = np.random.default_rng(0).integers(1900, 4000, IMG_SHAPE).astype(np.uint16) + raw = RawImage(stored, FILTER_ARRAY, np.array([2040.3, 2041.7, 2042.1, 2043.9]), 16383.0, request.param, + np.array([1.5, 2.5, 3.5, 4.5])) + monkeypatch.setattr(image_reading, '_read_raw_image', lambda filename: raw) + return ImageChannels('img.CR3') + + +def quantities_by_definition(img_channels, channel, row, col, sideband_pixels): + """The quantities of one LED, from their definitions on the arrays of the whole image.""" + data, mask, noise = img_channels.channel_data_mask_and_noise(channel) + box = np.s_[row - SEARCH_AREA_RADIUS:row + SEARCH_AREA_RADIUS, col - SEARCH_AREA_RADIUS:col + SEARCH_AREA_RADIUS] + values, in_channel = data[box], mask[box] + ring = np.ones(values.shape, dtype=bool) + ring[2:-2, 2:-2] = False + ring_pixels = values[ring & in_channel] + num_channel_pixels = np.count_nonzero(in_channel) + sideband_values = np.zeros(0) if sideband_pixels is None else data[sideband_pixels][mask[sideband_pixels]] + sideband_background = np.mean(sideband_values, dtype=np.float64) if sideband_values.size else np.nan + box_sum = float(np.sum(values)) + return {'sum_color_value': np.sum(values), 'mean_color_value': np.mean(values), + 'max_color_value': np.amax(values), + 'bgsub_sum_color_value': np.sum(values.astype(np.float64) - np.median(ring_pixels) * in_channel), + 'sidebg_sum_color_value': box_sum - num_channel_pixels * sideband_background, + 'sidebg_pedestal': num_channel_pixels * sideband_background, + 'bgsub2_sum_color_value': box_sum - num_channel_pixels * np.mean(ring_pixels, dtype=np.float64), + 'noise_sum_color_value': np.sqrt(np.sum(noise[box].astype(np.float64) ** 2))} + + +class TestStackedAnalysis: + @pytest.mark.parametrize('channel', [0, 1, 2]) + def test_quantities_are_those_of_their_definitions(self, img_channels, channel): + config = make_config() + led_data = generate_analysis_data('img.CR3', channel, SEARCH_AREAS, LINE_INDICES, config, fit_leds=False, + img_channels=img_channels) + sideband_pixels = sideband_pixels_of(img_channels, config) + assert [(data.led_id, data.led_array) for data in led_data] == [(0, 0), (1, 0), (2, 0), (3, 0), (4, 0), + (5, 1)] + for data, (led_id, row, col) in zip(led_data, SEARCH_AREAS): + expected = quantities_by_definition(img_channels, channel, row, col, sideband_pixels[led_id]) + for quantity in QUANTITIES: + assert getattr(data, quantity) == pytest.approx(expected[quantity], rel=1e-12, nan_ok=True), \ + f'{quantity} of LED {led_id}' + + def test_a_search_area_at_the_border_is_cut_short(self, img_channels): + # as data[row - R:row + R, col - R:col + R] cuts it; the search areas of the other shape are a stack of + # their own, so test_quantities_are_those_of_their_definitions covers both + analysis_pixels = _get_analysis_pixels(SEARCH_AREAS, SEARCH_AREAS[:, 0], sideband_pixels_of(img_channels, + make_config()), img_channels.shape, SEARCH_AREA_RADIUS) + stacks = {shape: positions.tolist() for positions, shape, _ in analysis_pixels.box_groups} + if img_channels.shape == IMG_SHAPE: + assert stacks == {(8, 8): [0, 1, 2, 3, 5], (6, 8): [4]} + else: + assert stacks == {(8, 8): [0, 1, 2, 3, 4], (8, 6): [5]} + + def test_the_led_without_side_bands_has_no_side_band_background(self, img_channels): + led_data = generate_analysis_data('img.CR3', 0, SEARCH_AREAS, LINE_INDICES, make_config(), fit_leds=False, + img_channels=img_channels) + assert np.isnan(led_data[5].sidebg_pedestal) + assert np.isfinite([data.sidebg_pedestal for data in led_data[:5]]).all() + + def test_only_the_analysed_pixels_are_calibrated_once_for_all_channels(self, img_channels, monkeypatch): + calibrated = [] + channel_data = image_reading._channel_data + monkeypatch.setattr(image_reading, '_channel_data', + lambda raw: calibrated.append(raw.data.size) or channel_data(raw)) + for channel in [0, 1, 2]: + generate_analysis_data('img.CR3', channel, SEARCH_AREAS, LINE_INDICES, make_config(), fit_leds=False, + img_channels=img_channels) + assert len(calibrated) == 1 + assert calibrated[0] < np.prod(IMG_SHAPE) + + def test_skipped_leds_are_not_analysed(self, img_channels): + config = make_config() + config['analyse_photo']['num_skip_leds'] = '1' + led_data = generate_analysis_data('img.CR3', 0, SEARCH_AREAS, LINE_INDICES, config, fit_leds=False, + img_channels=img_channels) + assert [data.led_id for data in led_data] == [0, 2, 4]