diff --git a/picamera_coverage.zip b/picamera_coverage.zip index 293a8ee6..fc4e7057 100644 Binary files a/picamera_coverage.zip and b/picamera_coverage.zip differ diff --git a/src/openflexure_microscope_server/background_detect.py b/src/openflexure_microscope_server/background_detect.py index d9e44cfc..9bd056dc 100644 --- a/src/openflexure_microscope_server/background_detect.py +++ b/src/openflexure_microscope_server/background_detect.py @@ -10,7 +10,6 @@ import cv2 import numpy as np from pydantic import BaseModel, Field, ConfigDict from pydantic.errors import PydanticUserError -from scipy.stats import norm from labthings_fastapi.thing_description import type_to_dataschema @@ -18,6 +17,14 @@ class MissingBackgroundDataError(RuntimeError): """An error raised if checking for sample without background data set.""" +class ChannelBlankError(RuntimeError): + """An error raised if a channel has no measured standard deviation. + + This is not physical and usually means the camera has not yet booted or changed + mode fully. + """ + + class BackgroundDetectorStatus(BaseModel): """The status information about a background detector instance needed for the GUI. @@ -136,7 +143,7 @@ class BackgroundDetectAlgorithm: "Each background detect algorithm must implement an image_is_sample method." ) - def set_background(self, image: np.ndarray) -> BaseModel: + def set_background(self, image: np.ndarray) -> None: """Use the input image to update the background data. Background data must be a Pydantic BaseModel. @@ -188,6 +195,10 @@ class ColourChannelDetectLUV(BackgroundDetectAlgorithm): background_data_model: BaseModel = ChannelDistributions settings_data_model: BaseModel = ColourChannelDetectSettings + # These are the same as those used for ChannelDeviationLUV. More detail is + # provided there. + min_stds = [0.5, 0.3, 0.5] + def background_mask(self, image: np.ndarray) -> np.ndarray: """Calculate a binary image, showing whether each pixel is background. @@ -196,11 +207,6 @@ class ColourChannelDetectLUV(BackgroundDetectAlgorithm): The image should be in LUV format, the output will be binary with the same shape in the first two dimensions. """ - if not self.background_data: - raise MissingBackgroundDataError( - "Background is not set: you need to calibrate background detection." - ) - # The ``[1:]`` selects only the U and V channels of the image. # Only U and V are used as brightness (L) often changes as # the height of the sample changes. @@ -225,6 +231,10 @@ class ColourChannelDetectLUV(BackgroundDetectAlgorithm): :returns: A value (between 0 and 100) is the percentage of the image that is sample. """ + if not self.background_data: + raise MissingBackgroundDataError( + "Background is not set: you need to calibrate background detection." + ) image_luv = cv2.cvtColor(image, cv2.COLOR_RGB2LUV) mask = self.background_mask(image_luv) @@ -250,15 +260,114 @@ class ColourChannelDetectLUV(BackgroundDetectAlgorithm): """Use the input image to update the background distributions.""" image_luv = cv2.cvtColor(image, cv2.COLOR_RGB2LUV) - ch1 = (image_luv.T[0]).flatten() - ch2 = (image_luv.T[1]).flatten() - ch3 = (image_luv.T[2]).flatten() + mu = np.mean(image_luv, axis=(0, 1)) + std = np.std(image_luv, axis=(0, 1)) - points = np.array([np.asarray(ch1), np.asarray(ch2), np.asarray(ch3)]).T - - # Get the mean and standard deviation of values in each channel - mu, std = np.apply_along_axis(norm.fit, 0, points) + if np.any(std == 0): + raise ChannelBlankError("Some LUV channels have no standard deviation.") + std = np.maximum(std, self.min_stds) self.background_data = ChannelDistributions( means=mu.tolist(), standard_deviations=std.tolist() ) + + +class ChannelDeviationLUV(BackgroundDetectAlgorithm): + """Compare the standard deviations of the LUV channels in a grid to background data. + + Using an LUV colour space, each image is divided into an 8x8 grid of images. + The standard deviation of each channel of each image is calculated and compared + to the median standard deviation for a grid of background images. + """ + + # Note we don't use the means in this algorithm but we use the same channel + # distributions model + background_data_model: BaseModel = ChannelDistributions + settings_data_model: BaseModel = ColourChannelDetectSettings + + # Empirically, 0.5 seems to be approximate the standard deviation for a good image + # in L and V. U appears to be about 60% of this value. U is about 65% of V when + # converting white-noise in RGB into LUV (note that we are using cv2's internal + # LUV colour space not converting to the CIELUV numbers) + min_stds = [0.5, 0.3, 0.5] + + def get_sample_coverage(self, image: np.ndarray) -> float: + """Return the percentage of the input image that is background. + + Evaluate whether it is foreground or background by comparing the standard + deviations of an 8x8 grid of sub-images to the median standard deviation + from a background image. + + :returns: A value (between 0 and 100) that is the percentage of the image that is + sample. + """ + if not self.background_data: + raise MissingBackgroundDataError( + "Background is not set: you need to calibrate background detection." + ) + image_luv = cv2.cvtColor(image, cv2.COLOR_RGB2LUV) + stds = _chunked_stds(image_luv, 8, 8) + + bg_stds = self.background_data.standard_deviations + l_cut = bg_stds[0] * self.settings.channel_tolerance + u_cut = bg_stds[1] * self.settings.channel_tolerance + v_cut = bg_stds[2] * self.settings.channel_tolerance + + populated_regions = ( + (stds[:, :, 0] > l_cut) | (stds[:, :, 1] > u_cut) | (stds[:, :, 2] > v_cut) + ) + + return float(100 * np.sum(populated_regions) / 64) + + def image_is_sample(self, image: np.ndarray) -> tuple[bool, str]: + """Label the current image as either background or sample. + + :returns: A tuple of the result (boolean), and explanation string. The + explanation string is formatted so it can be added into a sentence such as + ``An action was taken because the image is {message}.`` + """ + sample_coverage = self.get_sample_coverage(image) + + # Use bool otherwise get numpy variants of True and False. + is_sample = bool(sample_coverage > self.settings.min_sample_coverage) + message = f"{sample_coverage:0.1f}% sample" + if not is_sample: + message = "only " + message + return is_sample, message + + def set_background(self, image: np.ndarray) -> None: + """Use the input image to update the background distributions.""" + image_luv = cv2.cvtColor(image, cv2.COLOR_RGB2LUV) + mu = np.zeros(3) + c_stds = _chunked_stds(image_luv, 8, 8) + channel_blank = np.all(c_stds == 0, axis=(0, 1)) + if np.any(channel_blank): + raise ChannelBlankError("Some LUV channels have no standard devaition.") + + std = np.median(c_stds, axis=(0, 1)) + std = np.maximum(std, self.min_stds) + self.background_data = ChannelDistributions( + means=mu.tolist(), standard_deviations=std.tolist() + ) + + +def _chunked_stds(img: np.ndarray, n_rows: int = 8, n_cols: int = 8) -> np.ndarray: + """Split image into a grid and calculate std of each channel in each chunk. + + :param img: The image to analyse + :param n_rows: The number of rows in the grid + :param n_cols: The number of cols in the grid + :return: A numpy array of the grid of standard deviations. + """ + h, w = img.shape[:2] + row_height = h // n_rows + col_width = w // n_cols + out = np.zeros((n_rows, n_cols, 3)) + for i in range(n_rows): + for j in range(n_cols): + chunk = img[ + i * row_height : (i + 1) * row_height, + j * col_width : (j + 1) * col_width, + ] + out[i, j, :] = np.std(chunk, axis=(0, 1)) + return out diff --git a/src/openflexure_microscope_server/things/autofocus.py b/src/openflexure_microscope_server/things/autofocus.py index c8253f1d..c80b5308 100644 --- a/src/openflexure_microscope_server/things/autofocus.py +++ b/src/openflexure_microscope_server/things/autofocus.py @@ -28,6 +28,7 @@ from .stage import StageDependency as Stage LOGGER = logging.getLogger(__name__) MIN_TEST_IMAGE_COUNT = 3 MAX_TEST_IMAGE_COUNT = 9 +EXTRA_STACK_CAPTURES = 15 class NotStreamingError(RuntimeError): @@ -133,7 +134,7 @@ class StackParams(BaseModel): This is 15 images more then the minimum number that are captured. """ - return self.min_images_to_test + 15 + return self.min_images_to_test + EXTRA_STACK_CAPTURES def slice_to_save(self, sharpest_index: int) -> slice: """Return the slice of images to save given the index of the sharpest image.""" @@ -421,38 +422,38 @@ class AutofocusThing(lt.Thing): is close to focus, but not quite within ``dz/2``. It will attempt to autofocus up to 10 times. """ - repeat = True - attempts = 0 + attempt = 0 backlash = 200 with sharpness_monitor.run(): - while repeat and attempts < 10: + while attempt < 10: + attempt += 1 if start == "centre": stage.move_relative(x=0, y=0, z=-(backlash + dz / 2)) stage.move_relative(x=0, y=0, z=backlash) + # Always start centrally for future runs + start = "centre" + focus_data_index, _ = sharpness_monitor.focus_rel( dz, block_cancellation=True ) _times, heights, sizes = sharpness_monitor.move_data(focus_data_index) peak_height = heights[np.argmax(sizes)] - height_min = np.min(heights) - height_max = np.max(heights) + target_min = np.min(heights) + dz / 5 + target_max = np.max(heights) - dz / 5 - if ( - peak_height - height_min < dz / 5 - or height_max - peak_height < dz / 5 - ): - attempts += 1 - start = "centre" - stage.move_absolute(z=peak_height - backlash) - stage.move_absolute(z=peak_height) - else: - repeat = False - stage.move_relative(x=0, y=0, z=-(dz + backlash)) - stage.move_absolute(z=peak_height) - return heights.tolist(), sizes.tolist() + # move to the peak + stage.move_absolute(z=peak_height - backlash) + stage.move_absolute(z=peak_height) + + if target_min < peak_height < target_max: + # If it is within the target range then return + return heights.tolist(), sizes.tolist() + raise NoFocusFoundError( + "Looping autofocus couldn't converge on a focus location." + ) stack_images_to_save = lt.ThingSetting( initial_value=1, @@ -566,7 +567,9 @@ class AutofocusThing(lt.Thing): stage: Stage, sharpness_monitor: SharpnessMonitorDep, stack_parameters: StackParams, - ) -> tuple[bool, int]: + save_on_failure: bool = False, + check_turning_points: bool = True, + ) -> tuple[bool, Optional[int]]: """Run a smart stack. A smart stack captures images offset in z, testing whether the sharpest image @@ -582,17 +585,22 @@ class AutofocusThing(lt.Thing): supplied by LabThings dependency injection :param stack_parameters: A StackParams object containing the required parameters to run a stack. + :param save_on_failure: Whether to save an image even if no focus was found. + :param check_turning_points: Whether to check the number of turning points in + the sharpnesses of the images in the stack is exactly 1. (May fail with + thick samples) :returns: A tuple containing: * A boolean, True if stack was successfully * The z position of the sharpest image """ - tries = 0 - # Loop until a stack is successful - while tries < stack_parameters.max_attempts: + attempt = 0 + while True: + attempt += 1 success, captures, sharpest_id = self.z_stack( stack_parameters=stack_parameters, + check_turning_points=check_turning_points, cam=cam, stage=stage, ) @@ -600,28 +608,33 @@ class AutofocusThing(lt.Thing): if success: break + if attempt >= stack_parameters.max_attempts: + break + # The z position of the first images in the previous attempt. initial_z_pos = captures[0].position["z"] # If a stack is not successful, move to the start and autofocus - self.reset_stack( - initial_z_pos, - stack_parameters.autofocus_dz, - stage, - sharpness_monitor, - ) + try: + self.reset_stack( + initial_z_pos, + stack_parameters.autofocus_dz, + stage, + sharpness_monitor, + ) + except NoFocusFoundError: + break # Save stack_parameters.image_to_save images centred on the sharpest capture. # If the smart_stack failed the exact number of images saved may not be # stack_parameters.image_to_save - self.save_stack( - sharpest_id=sharpest_id, - captures=captures, - stack_parameters=stack_parameters, - cam=cam, - ) + if success or save_on_failure: + self.save_stack( + sharpest_id=sharpest_id, + captures=captures, + stack_parameters=stack_parameters, + cam=cam, + ) - # Return whether or not the smart stack was successful, and the z position of - # the sharpest image, for path planning and tracking return success, _get_capture_by_id(captures, sharpest_id).position["z"] def reset_stack( @@ -681,9 +694,10 @@ class AutofocusThing(lt.Thing): def z_stack( self, stack_parameters: StackParams, + check_turning_points: bool, cam: CameraClient, stage: Stage, - ) -> tuple[bool, list[CaptureInfo], Optional[int]]: + ) -> tuple[bool, list[CaptureInfo], int]: """Capture a series of images checking that sharpest image central. The images are separated in z offset by stack_parameters.stack_dz, as they @@ -692,6 +706,9 @@ class AutofocusThing(lt.Thing): completes. :param stack_parameters: a StackParams object holding stack parameters + :param check_turning_points: Whether to check the number of turning points in + the sharpnesses of the images in the stack is exactly 1. (May fail with + thick samples) :param cam: Camera Dependency to be passed through from the calling action :param stage: Stage Dependency to be passed through from the calling action @@ -699,7 +716,7 @@ class AutofocusThing(lt.Thing): * the stack result (True for successful stack, False for failed stack), * a list of CaptureInfo objects, - * the buffer_id of the sharpest image (or None if the stack failed). + * the buffer_id of the sharpest image """ # Move down by the height of the z stack, plus an overshoot # Better to start too low and take too many images than too high and need to refocus @@ -718,7 +735,7 @@ class AutofocusThing(lt.Thing): ims_to_check = slice(-stack_parameters.min_images_to_test, None) # If the sharpest image isn't found within the maximum number of images - # end the loop and return "restart" + # end the loop and return False indicating the stack failed. while len(captures) < stack_parameters.max_images_to_test: time.sleep(stack_parameters.settling_time) @@ -733,16 +750,18 @@ class AutofocusThing(lt.Thing): # If the number of images is enough to test, test them if len(captures) >= stack_parameters.min_images_to_test: - result, capture_id = self.check_stack_result(captures[ims_to_check]) + result, capture_id = self.check_stack_result( + captures[ims_to_check], check_turning_points=check_turning_points + ) if result == "success": return True, captures, capture_id if result == "restart": - return False, captures, None + return False, captures, capture_id # If reached here the result was "continue" stage.move_relative(z=stack_parameters.stack_dz) - return False, captures, None + return False, captures, capture_id def capture_stack_image( self, @@ -774,12 +793,15 @@ class AutofocusThing(lt.Thing): # unlikely to improve readability. This function is basically a complex switch # statement, having an explicit return after each option is clear. def check_stack_result( # noqa: PLR0911 - self, captures: list[CaptureInfo] + self, captures: list[CaptureInfo], check_turning_points: bool ) -> tuple[Literal["success", "continue", "restart"], int]: """Check if the sharpest image in a list of captures is central enough. :param captures: a list of the capture objects to for testing if the sharpness has converged in the centre + :param check_turning_points: Whether to check the number of turning points in + the sharpnesses of the images in the stack is exactly 1. (May fail with + thick samples) :returns: A tuple with two values: @@ -793,28 +815,91 @@ class AutofocusThing(lt.Thing): * capture_id - the buffer id of the sharpest image """ - sharpest_index = np.argmax([capture.sharpness for capture in captures]) + sharpnesses = np.array([capture.sharpness for capture in captures]) + sharpest_index = np.argmax(sharpnesses) # The buffer id of the sharpest image capture_id = captures[sharpest_index].buffer_id - sharpness_length = len(captures) + n_imgs = len(captures) # If only testing one image, then by definition the sharpest is central - if sharpness_length == 1: + if n_imgs == 1: return "success", capture_id # If testing three images, test if the centre is the sharpest - if sharpness_length == 3: + if n_imgs == 3: if sharpest_index == 1: return "success", capture_id if sharpest_index == 0: return "restart", capture_id return "continue", capture_id - # For larger stacks, test if the best image is not within two of the edge of the stack - # ie for a stack of 7 images, best image must be between 3rd and and 5th - exclusion_range = 2 - - if sharpest_index < exclusion_range: - return "restart", capture_id - if sharpest_index >= sharpness_length - exclusion_range: + try: + turning = _get_peak_turning_point(sharpnesses) + except NotAPeakError: return "continue", capture_id + + # For larger stacks, test if the best image is not within two of the edge of + # the stack ie for a stack of 7 images, best image must be between 3rd and 5th + edge_size = 2 + + # Check both the peak from fitting and the sharpest image are not at the + # edge of the stack. + if turning < edge_size - 0.5 or sharpest_index < edge_size: + return "restart", capture_id + if turning > n_imgs - edge_size - 0.5 or sharpest_index >= n_imgs - edge_size: + return "continue", capture_id + + if check_turning_points: + turning_points = _count_turning_points(sharpnesses) + if turning_points > 1: + return "continue", capture_id + return "success", capture_id + + +class NotAPeakError(RuntimeError): + """The data to fit isn't a peak.""" + + +class NoFocusFoundError(RuntimeError): + """No focus found during looping Autofocus.""" + + +def _get_peak_turning_point(sharpnesses: np.ndarray) -> float: + """Get the turning point for a sharpnesses in a z-stack. + + :param sharpnesses: A numpy array of sharpnesses + :return: The x value of the turning point where x-axis is 0 to N-1 for the N + sharpness values + + :raise NotAPeakError: If the fit doesn't have a maximum within 95% confidence. + """ + # Fit the peak + x = range(len(sharpnesses)) + coeffs, cov = np.polyfit(x, sharpnesses, deg=2, cov=True) + # a in the equation a*x^2 + bx + c + a = float(coeffs[0]) + fit_func = np.poly1d(coeffs) + + # find the peaks's x-position + turning = float(fit_func.deriv().roots[0]) + # sigma_a the standard error of a + sigma_a = float(np.sqrt(cov[0, 0])) + # Estimate the upper 95% confidence bound for the x^2 term + ci_high = a + 2 * sigma_a + # If the high 95% confindence interval of the peak is not negative then the + # fit isn't sure if this is a peak or a U shape, so we continue. + # -1e-9 is used to guard against weird effects for flat and straight data + if ci_high > -1e-9: + raise NotAPeakError("Not a peak to within 95% confidence.") + return turning + + +def _count_turning_points(sharpnesses: np.ndarray) -> int: + """Count the number of turing points, after rejecting those from noise.""" + # Also take the difference of neighbouring points to check for turning points + d_sharpnesses = sharpnesses[1:] - sharpnesses[:-1] + # Filter out any points that are not prominent + prominent = abs(d_sharpnesses) / np.mean(abs(d_sharpnesses)) > 0.5 + d_sharpnesses = d_sharpnesses[prominent] + # count the sign changes + return int(np.sum(d_sharpnesses[1:] * d_sharpnesses[:-1] < 0)) diff --git a/src/openflexure_microscope_server/things/camera/__init__.py b/src/openflexure_microscope_server/things/camera/__init__.py index fc11ca1f..1bb3bcef 100644 --- a/src/openflexure_microscope_server/things/camera/__init__.py +++ b/src/openflexure_microscope_server/things/camera/__init__.py @@ -28,6 +28,7 @@ from labthings_fastapi.types.numpy import NDArray from openflexure_microscope_server.ui import ActionButton, PropertyControl from openflexure_microscope_server.background_detect import ( ColourChannelDetectLUV, + ChannelDeviationLUV, BackgroundDetectAlgorithm, BackgroundDetectorStatus, ) @@ -182,8 +183,11 @@ class BaseCamera(lt.Thing): dictionary in this function. Configuration will be added at a later date. """ super().__init__() - self.background_detectors = {"Colour Channels (LUV)": ColourChannelDetectLUV()} - self._detector_name = "Colour Channels (LUV)" + self.background_detectors = { + "Colour Channels (LUV)": ColourChannelDetectLUV(), + "Channel Deviations (LUV)": ChannelDeviationLUV(), + } + self._detector_name = "Channel Deviations (LUV)" def __enter__(self) -> None: """Open hardware connection when the Thing context manager is opened.""" @@ -637,7 +641,7 @@ class BaseCamera(lt.Thing): def detector_name(self, name: str) -> None: """Validate and set detector_name.""" if name not in self.background_detectors: - raise ValueError(f"{name} is not a valid background detector name") + LOGGER.warning(f"{name} is not a valid background detector name.") self._detector_name = name @property diff --git a/src/openflexure_microscope_server/things/camera/picamera.py b/src/openflexure_microscope_server/things/camera/picamera.py index ff7c8867..7c664dbc 100644 --- a/src/openflexure_microscope_server/things/camera/picamera.py +++ b/src/openflexure_microscope_server/things/camera/picamera.py @@ -42,6 +42,7 @@ from openflexure_microscope_server.ui import ( action_button_for, property_control_for, ) +from openflexure_microscope_server.background_detect import ChannelBlankError from . import picamera_recalibrate_utils as recalibrate_utils from . import picamera_tuning_file_utils as tf_utils @@ -767,7 +768,6 @@ class StreamingPiCamera2(BaseCamera): * ``auto_expose_from_minimum`` * ``set_static_green_equalisation`` to set geq offset to max * ``calibrate_lens_shading`` (also sets colour gains for white balance) - * ``set_background`` """ self.flat_lens_shading() @@ -775,8 +775,16 @@ class StreamingPiCamera2(BaseCamera): self.set_static_green_equalisation() self.set_ce_enable_to_off() self.calibrate_lens_shading() - time.sleep(0.5) - self.set_background(portal) + for _i in range(3): + try: + time.sleep(self._sensor_info.long_pause) + self.set_background(portal) + # Return if background is set + return + except ChannelBlankError: + # If channel is blank, sleep a second and try again. + pass + raise RuntimeError("Couldn't set background") @lt.thing_property def primary_calibration_actions(self) -> list[ActionButton]: diff --git a/src/openflexure_microscope_server/things/smart_scan.py b/src/openflexure_microscope_server/things/smart_scan.py index c9454b11..0d35cc2f 100644 --- a/src/openflexure_microscope_server/things/smart_scan.py +++ b/src/openflexure_microscope_server/things/smart_scan.py @@ -467,13 +467,17 @@ class SmartScanThing(lt.Thing): continue focused, focused_height = self._autofocus.run_smart_stack( - stack_parameters=self._stack_params + stack_parameters=self._stack_params, + save_on_failure=not self._scan_data.skip_background, ) current_pos_xyz = (new_pos_xyz[0], new_pos_xyz[1], focused_height) + # An image was captured if we are focussed or we are not skipping background. + imaged = focused or not self._scan_data.skip_background + route_planner.mark_location_visited( - current_pos_xyz, imaged=True, focused=focused + current_pos_xyz, imaged=imaged, focused=focused ) # increment capture counter as thread has completed diff --git a/tests/test_autofocus.py b/tests/test_autofocus.py new file mode 100644 index 00000000..b68dfb1a --- /dev/null +++ b/tests/test_autofocus.py @@ -0,0 +1,97 @@ +"""Tests for the autofoucs logic. + +This doesn't check the behaviour of the JPEG shaprness monitor. +""" + +import pytest +import numpy as np + +from openflexure_microscope_server.things.autofocus import ( + AutofocusThing, + NoFocusFoundError, +) + + +def fake_sharpness_data( + dz: int, start_z: int, max_loc: int, length: int = 41 +) -> tuple[list[float], np.ndarray, np.ndarray]: + """Create some fake data for the shapeness. + + The highest returned sharpness is closest to max_loc + """ + # Some fake timestamps + times = [i / 10 + 100000 for i in range(length)] + img_dz = dz / (length - 1) + heights = [round(start_z + i * img_dz) for i in range(length)] + # Sharpnesses fall off linearly in this model. + sharpnesses = [10 * dz - abs(max_loc - h) for h in heights] + return times, np.array(heights), np.array(sharpnesses) + + +@pytest.mark.parametrize( + ("start_z", "max_loc", "centre", "attempts_expected", "passes"), + [ + # To complete, the max must be in the central 1200, so -600 to 600 when looping + # from -1000 to 1000 + (0, 550, True, 1, True), # Found in loop1 from -1000 to 1000 + (0, 650, True, 2, True), # Just outside the limit in loop1 + (0, 1300, True, 2, True), # Found in loop2 from 0 to 2000 + (0, 1300, False, 1, True), # Found in loop1 from 0 to 2000 (as start="base") + (0, -1300, True, 2, True), # Found in loop2 from 0 to 2000 + (0, 7300, True, 8, True), # Found in loop8 from 6000 to 8000 + (0, 9300, True, 10, True), # Found in loop10 from 8000 to 10000 + (0, 9900, True, 10, False), # Still not central in loop 10, doesn't pass + ], +) +def test_looping_autofocus(start_z, max_loc, centre, attempts_expected, passes, mocker): + """Test the high level looping autofocus algorithm.""" + dz = 2000 + # Make a mock stage where move_absolute abs and relative updates the position counter. + stage = mocker.Mock() + stage.position = {"x": 0, "y": 0, "z": start_z} + + def set_pos(**kwargs: int) -> None: + """Move absolute should update position. So make a side effect for the mock.""" + for axis, value in kwargs.items(): + stage.position[axis] = value + + def adjust_pos(**kwargs: int) -> None: + """Move relative should update position. So make a side effect for the mock.""" + for axis, value in kwargs.items(): + stage.position[axis] += value + + stage.move_absolute.side_effect = set_pos + stage.move_relative.side_effect = adjust_pos + + # Make a mock sharpness monitor that can generate sharpness data. + sharpness_monitor = mocker.MagicMock() + sharpness_monitor.focus_rel.return_value = (0, 0) + + def return_sharpness(*_args) -> tuple[list[float], np.ndarray, np.ndarray]: + """Generate sharpnesses based on parameterised input, and mock stage position.""" + return fake_sharpness_data( + dz=dz, + start_z=stage.position["z"], + max_loc=max_loc, + ) + + sharpness_monitor.move_data.side_effect = return_sharpness + + autofocus_thing = AutofocusThing() + if passes: + autofocus_thing.looping_autofocus( + stage=stage, + sharpness_monitor=sharpness_monitor, + dz=dz, + start="centre" if centre else "base", + ) + else: + with pytest.raises(NoFocusFoundError): + autofocus_thing.looping_autofocus( + stage=stage, + sharpness_monitor=sharpness_monitor, + dz=dz, + start="centre" if centre else "base", + ) + assert sharpness_monitor.focus_rel.call_count == attempts_expected + assert abs(max_loc - stage.position["z"]) < dz / 40 diff --git a/tests/test_background_detectors.py b/tests/test_background_detectors.py index 956a96e9..37df95ae 100644 --- a/tests/test_background_detectors.py +++ b/tests/test_background_detectors.py @@ -15,6 +15,9 @@ from openflexure_microscope_server.background_detect import ( ChannelDistributions, ColourChannelDetectSettings, ColourChannelDetectLUV, + _chunked_stds, + ChannelDeviationLUV, + ChannelBlankError, ) RNG = np.random.default_rng() @@ -181,3 +184,140 @@ def test_colour_channel_luv_load_bad_data(): cc_luv.settings = WrongModel() with pytest.raises(TypeError): cc_luv.background_data = WrongModel() + + +def create_patchwork_image(magnitude=3, blank_channels=None): + """Create a patchwork image, with known stds per chunk. + + :return: The image, and the (8,8,3) array of stds + """ + if blank_channels is None: + blank_channels = [] + n_rows = 8 + n_cols = 8 + chunk_h = 10 + chunk_w = 10 + + # Precompute chunk stds and construct the final image + expected_stds = np.zeros((n_rows, n_cols, 3), dtype=float) + img = np.zeros((n_rows * chunk_h, n_cols * chunk_w, 3), dtype=np.uint8) + + for i in range(n_rows): + for j in range(n_cols): + for channel in range(3): + if channel in blank_channels: + chunk = np.zeros((chunk_h, chunk_w), dtype=np.uint8) + else: + chunk = 9.8 * np.ones((chunk_h, chunk_w)) + chunk += RNG.normal(scale=magnitude, size=(chunk_h, chunk_w)) + chunk = chunk.astype(np.uint8) + + # Store its std + expected_stds[i, j, channel] = np.std(chunk) + + # Insert chunk into the final large image + y0, y1 = i * chunk_h, (i + 1) * chunk_h + x0, x1 = j * chunk_w, (j + 1) * chunk_w + img[y0:y1, x0:x1, channel] = chunk + return img, expected_stds + + +def test_chunked_stds_with_precomputed_chunk_stds(): + """Test _chunked_stds returns the answer calculated when making a patchwork image.""" + img, expected_stds = create_patchwork_image() + + # Run the function under test + result = _chunked_stds(img, n_rows=8, n_cols=8) + + # Compare + np.testing.assert_allclose(result, expected_stds, rtol=1e-6, atol=1e-12) + + +def test_channel_deviation_luv_set_background(mocker): + """Test set_background takes the median of each channel, and errors for blank channels.""" + cd_luv = ChannelDeviationLUV() + + # Patch RGB to LUV so or we don't know what the STDs should be + mocker.patch("cv2.cvtColor", side_effect=lambda img, _method: img) + + for magnitude in [0.2, 1, 10]: + # Create an image and set it as background + img, expected_stds = create_patchwork_image(magnitude=magnitude) + cd_luv.set_background(img) + + # Do a somewhat verbose checking for clarity + for channel in range(3): + # Saved std + channel_std = cd_luv.background_data.standard_deviations[channel] + # Expected median + channel_median = np.median(expected_stds[:, :, channel]) + # If the median is above the minimum allowed then it should be returned + if channel_median > cd_luv.min_stds[channel]: + assert np.isclose(channel_std, channel_median, rtol=1e-6, atol=1e-12) + else: + # If not the minimum is returned. + assert channel_std == cd_luv.min_stds[channel] + + # Also check that if channels are blank then an error is thrown + img, _expected_stds = create_patchwork_image(magnitude=1, blank_channels=[1, 2]) + with pytest.raises(ChannelBlankError): + cd_luv.set_background(img) + + +def test_channel_deviation_luv_image_is_sample(background_image, mocker): + """Check image_is_sample reports the result from get_sample_coverage.""" + cd_luv = ChannelDeviationLUV() + + # No background data so it is not ready and will error if image_is_sample is called. + assert not cd_luv.status.ready + with pytest.raises(MissingBackgroundDataError): + cd_luv.image_is_sample(background_image) + + cd_luv.settings.min_sample_coverage = 20 + cd_luv.get_sample_coverage = mocker.Mock(return_value=10) + is_sample, message = cd_luv.image_is_sample(background_image) + assert not is_sample + assert message == r"only 10.0% sample" + + # Reduce the min coverage + cd_luv.settings.min_sample_coverage = 9 + + is_sample, message = cd_luv.image_is_sample(background_image) + assert is_sample + assert message == r"10.0% sample" + + +def test_channel_deviation_luv_get_sample_coverage(background_image, mocker): + """Check _get_sample_coverage returns the values expected.""" + cd_luv = ChannelDeviationLUV() + + # Create fake chunked STD data where each channel is the numbers 0 -> 31.5 in 0.5 + # steps + grid = np.arange(0, 32, 0.5).reshape(8, 8) + fake_stds = np.stack([grid, grid, grid], axis=-1) + mocker.patch( + "openflexure_microscope_server.background_detect._chunked_stds", + return_value=fake_stds, + ) + + # Create fake background + cd_luv.background_data = ChannelDistributions( + means=[0, 0, 0], standard_deviations=[1.1, 1.1, 1.1] + ) + # Get sample coverage with channel tolerance of 7. Checking each channel for the + # numbers below 7.7. There are 16 out of 64. So 75% should be sample + cd_luv.settings.channel_tolerance = 7 + assert cd_luv.get_sample_coverage(background_image) == 75 + # This is unchanged if two channels have larger background values. + cd_luv.background_data = ChannelDistributions( + means=[0, 0, 0], standard_deviations=[1.6, 1.6, 1.1] + ) + assert cd_luv.get_sample_coverage(background_image) == 75 + # But coverage increases if any channels has a lower background value. + cd_luv.background_data = ChannelDistributions( + means=[0, 0, 0], standard_deviations=[1.6, 0.6, 1.1] + ) + assert cd_luv.get_sample_coverage(background_image) == 85.9375 + # Returns to 75% if that channel is empty + fake_stds[:, :, 1] = 0 + assert cd_luv.get_sample_coverage(background_image) == 75 diff --git a/tests/test_stack.py b/tests/test_stack.py index 747e2eb1..57e493db 100644 --- a/tests/test_stack.py +++ b/tests/test_stack.py @@ -1,7 +1,4 @@ -"""Tests for the smart/fast stacking. - -Currently these tests don't test the Thing itself, just surrounding functionality -""" +"""Tests for the smart and fast stacking.""" from typing import Optional import tempfile @@ -23,6 +20,10 @@ from openflexure_microscope_server.things.autofocus import ( MAX_TEST_IMAGE_COUNT, _get_capture_by_id, _get_capture_index_by_id, + EXTRA_STACK_CAPTURES, + NotAPeakError, + _get_peak_turning_point, + _count_turning_points, ) from openflexure_microscope_server.scan_directories import IMAGE_REGEX @@ -352,3 +353,338 @@ def test_coercing_stack_save_ims( assert stack_params.images_to_save == coerced_save_ims # Check that the setting in the Thing was updated to the coerced value assert stack_params.images_to_save == autofocus_thing.stack_images_to_save + + +@pytest.mark.parametrize("pass_on", [1, 2, 3, 4]) +def test_run_smart_stack(pass_on, autofocus_thing, mocker): + """Test Running smart stack with the stack passing on different attempts.""" + cam = mocker.Mock() + stage = mocker.Mock() + sharpness_monitor = mocker.MagicMock() + stack_params = autofocus_thing.create_stack_params( + autofocus_dz=2000, + images_dir="/this/is/fake", + save_resolution=(1640, 1232), + logger=LOGGER, + ) + assert stack_params.max_attempts == 3 + + # Set up returns from z-stack + fake_captures = [ + CaptureInfo( + buffer_id="first", position={"x": 0, "y": 0, "z": -99}, sharpness=123 + ), + CaptureInfo( + buffer_id="pick_me", position={"x": 0, "y": 0, "z": 555}, sharpness=456 + ), + CaptureInfo( + buffer_id="last", position={"x": 0, "y": 0, "z": 999}, sharpness=123 + ), + ] + + successful_return = (True, fake_captures, "pick_me") + failed_return = (False, fake_captures, "pick_me") + return_list = [failed_return] * (pass_on - 1) + [successful_return] + + # Mock z_stack and looping_autofocus + autofocus_thing.z_stack = mocker.Mock(side_effect=return_list) + autofocus_thing.looping_autofocus = mocker.Mock() + + # Run it + success, final_z = autofocus_thing.run_smart_stack( + cam=cam, + stage=stage, + sharpness_monitor=sharpness_monitor, + stack_parameters=stack_params, + save_on_failure=False, + check_turning_points=True, + ) + + # Only passes if the attempt it passes on is less than max attempts + assert success == (pass_on <= stack_params.max_attempts) + # Final z is the one from the id returned by the stack "pick_me" + assert final_z == 555 + + # z_stack should run up until the time it passes. Running no more than max_attempts + n_stacks = min(pass_on, stack_params.max_attempts) + assert autofocus_thing.z_stack.call_count == n_stacks + # Move absolute should be 1 less time that the number of times z_stack_run + assert stage.move_absolute.call_count == n_stacks - 1 + # As should looping autofocus + assert autofocus_thing.looping_autofocus.call_count == n_stacks - 1 + + # Check rest stack is moving to the first image in the stack. + if n_stacks > 1: + assert stage.move_absolute.call_args.kwargs["z"] == -99 + + # Mock called to save image + assert cam.save_from_memory.call_count == (1 if success else 0) + + +def setup_and_run_z_stack(check_returns, check_turning_points, autofocus_thing, mocker): + """Set up a z_stack, run it, and return the result. + + :param check_returns: The return values from check_stack_result. Note that if this + is a list, it will be set as a side effect (and should be a list of tuples of + results). If it a tuple (or anything else), it is set as a return value. + """ + stack_params = autofocus_thing.create_stack_params( + autofocus_dz=2000, + images_dir="/this/is/fake", + save_resolution=(1640, 1232), + logger=LOGGER, + ) + stack_params.settling_time = 0 # Don't settle or tests take forever. + + stage = mocker.Mock() + cam = mocker.Mock() + autofocus_thing.capture_stack_image = mocker.Mock() + if isinstance(check_returns, list): + autofocus_thing.check_stack_result = mocker.Mock(side_effect=check_returns) + else: + autofocus_thing.check_stack_result = mocker.Mock(return_value=check_returns) + return autofocus_thing.z_stack( + stack_parameters=stack_params, + check_turning_points=check_turning_points, + cam=cam, + stage=stage, + ) + + +def test_z_stack_turning_toggle_passed(autofocus_thing, mocker): + """Check that the toggling of turning points is passed to the check.""" + check_returns = ("success", "mock_id") + for check_turning in [True, False]: + setup_and_run_z_stack(check_returns, check_turning, autofocus_thing, mocker) + check_kwargs = autofocus_thing.check_stack_result.call_args.kwargs + assert check_kwargs["check_turning_points"] == check_turning + + +def test_z_stack_returns_on_success_and_restart(autofocus_thing, mocker): + """Check that if the check returns success or restart then the stack exits with correct return value.""" + for result in ["success", "restart"]: + check_returns = (result, "mock_id") + ret = setup_and_run_z_stack(check_returns, True, autofocus_thing, mocker) + assert autofocus_thing.check_stack_result.call_count == 1 + # Check the number of images taken is exactly the call count. + ims_taken = autofocus_thing.capture_stack_image.call_count + assert ims_taken == autofocus_thing.stack_min_images_to_test + # And the result is as expected. + assert ret[0] == (result == "success") + + +def test_z_stack_exits_if_focus_never_found(autofocus_thing, mocker): + """Check that if the check returns continue the stack exits eventually with a failure.""" + check_returns = ("continue", "mock_id") + ret = setup_and_run_z_stack(check_returns, True, autofocus_thing, mocker) + + assert autofocus_thing.check_stack_result.call_count == EXTRA_STACK_CAPTURES + 1 + # Check the number of images taken is the maximum possible, set by the min images to + # test and the number of extra images that can be taken + ims_taken = autofocus_thing.capture_stack_image.call_count + max_ims = autofocus_thing.stack_min_images_to_test + EXTRA_STACK_CAPTURES + assert ims_taken == max_ims + # And the result is as expected. + assert not ret[0] + + +def test_z_stack_return(autofocus_thing, mocker): + """Check z-stack returns as expected for more complex cases the fixed results above.""" + for i in range(2, EXTRA_STACK_CAPTURES): + check_returns = [ + ("restart" if j == i - 1 else "continue", f"id_{j}") for j in range(i) + ] + ret = setup_and_run_z_stack(check_returns, True, autofocus_thing, mocker) + # Calculate images taken + images_taken = autofocus_thing.stack_min_images_to_test + i - 1 + assert autofocus_thing.capture_stack_image.call_count == images_taken + # Check it reports a failure + assert not ret[0] + + # Repeat ending with a success rather than a failure + check_returns = [ + ("success" if j == i - 1 else "continue", f"id_{j}") for j in range(i) + ] + ret = setup_and_run_z_stack(check_returns, True, autofocus_thing, mocker) + # Calculate images taken + assert autofocus_thing.capture_stack_image.call_count == images_taken + # Check it reports a success + assert ret[0] + + +def test_capture_stack_image(autofocus_thing, mocker): + """Check that capture stack image calls the expected functions and returns the expected data.""" + stage = mocker.Mock() + stage.position = {"x": 123, "y": 456, "z": 789} + cam = mocker.Mock() + cam.capture_to_memory.return_value = "fake_buffer_id" + cam.grab_jpeg_size.return_value = 54321 + buffer_max = 11 + + info = autofocus_thing.capture_stack_image( + cam=cam, stage=stage, buffer_max=buffer_max + ) + assert cam.capture_to_memory.call_count == 1 + assert cam.grab_jpeg_size.call_count == 1 + assert info.buffer_id == "fake_buffer_id" + assert info.position == {"x": 123, "y": 456, "z": 789} + assert info.sharpness == 54321 + + +def mock_capture(buffer_id: int, sharpness: int) -> CaptureInfo: + """Create a CaptureInfo instance with a dummy position.""" + return CaptureInfo( + buffer_id=buffer_id, + position={"x": 0, "y": 0, "z": buffer_id}, + sharpness=sharpness, + ) + + +def test_check_stack_single_image_returns_success(autofocus_thing): + """A single image is always successful.""" + captures = [mock_capture("mock-id", 10)] + result, cap_id = autofocus_thing.check_stack_result( + captures, check_turning_points=False + ) + assert result == "success" + assert cap_id == "mock-id" + + +@pytest.mark.parametrize( + ("sharpnesses", "expected"), + [ + ([5, 10, 3], "success"), + ([10, 4, 2], "restart"), + ([1, 2, 10], "continue"), + ], +) +def test_check_stack_three_image_logic(sharpnesses, expected, autofocus_thing): + """For 3 images, success is the highest one is central.""" + captures = [mock_capture(i, s) for i, s in enumerate(sharpnesses)] + result, _ = autofocus_thing.check_stack_result(captures, check_turning_points=False) + assert result == expected + + +def _run_check_stack_with_good_peak(autofocus_thing, count_turnings=False): + """Run check stack on a good peak that should pass, and return the result. + + This can be used to check how other mocked results of subfunctions affects the + result. + """ + # Create an obvious peak that would normally pass. + sharpnesses = [1, 2, 4, 7, 12, 7, 4, 2, 1] + captures = [mock_capture(i, s) for i, s in enumerate(sharpnesses)] + + result, cap_id = autofocus_thing.check_stack_result( + captures, check_turning_points=count_turnings + ) + # Nothing a mocked function does should change which is the sharpest image. + assert cap_id == 4 + return result + + +def test_check_stack_continues_if_no_tuning_point(autofocus_thing, mocker): + """Check that continue is returned if no turning point is found.""" + # Mock to simulate not finding a peak + mocker.patch( + "openflexure_microscope_server.things.autofocus._get_peak_turning_point", + side_effect=NotAPeakError, + ) + result = _run_check_stack_with_good_peak(autofocus_thing) + # Check that the NotAPeakError causes it to continue instead. + assert result == "continue" + + +@pytest.mark.parametrize( + ("location", "expected"), + [ + (-10, "restart"), # Restart if lower than 1.5 (halfway between im 2 and 3) + (-1, "restart"), + (0, "restart"), + (1, "restart"), + (1.49, "restart"), + (1.5, "success"), # Success up to 6.5 (as we have 9 images, final index is 8) + (2.5, "success"), + (4.5, "success"), + (6.5, "success"), + (6.51, "continue"), # Continue if thrung point is after 6.5 + (7, "continue"), + (8.1, "continue"), + (123, "continue"), + ], +) +def test_check_stack_affected_by_turning_point_location( + location, expected, autofocus_thing, mocker +): + """Check that the turning point location affects the return as expected.""" + # Mock to give the turning point location specified + mocker.patch( + "openflexure_microscope_server.things.autofocus._get_peak_turning_point", + return_value=location, + ) + result = _run_check_stack_with_good_peak(autofocus_thing) + assert result == expected + + +def test_check_stack_affected_by_number_of_turning_points(autofocus_thing, mocker): + """Check that the turning point location affects the return as expected.""" + # Set the turning point to the centre + mocker.patch( + "openflexure_microscope_server.things.autofocus._get_peak_turning_point", + return_value=5, + ) + mocker.patch( + "openflexure_microscope_server.things.autofocus._count_turning_points", + return_value=1, + ) + result = _run_check_stack_with_good_peak(autofocus_thing, count_turnings=True) + # Successful with 1 peak + assert result == "success" + + # Change return to be 2 peaks + mocker.patch( + "openflexure_microscope_server.things.autofocus._count_turning_points", + return_value=2, + ) + result = _run_check_stack_with_good_peak(autofocus_thing, count_turnings=True) + # Continue with 2 peaks + assert result == "continue" + # Unless this check is turned off + result = _run_check_stack_with_good_peak(autofocus_thing, count_turnings=False) + assert result == "success" + + +def test_get_peak_turning_point(): + """Check that the peak fitting returns expected value (or error).""" + with pytest.raises(NotAPeakError): + _get_peak_turning_point(np.ones(9)) + + linear = np.arange(9) + u_shape = 2 * (linear - 4) ** 2 + 17 + peak = -2 * (linear - 4) ** 2 + 55 + + with pytest.raises(NotAPeakError): + _get_peak_turning_point(linear) + + with pytest.raises(NotAPeakError): + _get_peak_turning_point(u_shape) + + # Should be 4 to within a fitting error + assert abs(_get_peak_turning_point(peak) - 4) < 1e-7 + + +def test_count_turning_points(): + """Check the turing point count works as expected.""" + linear = np.arange(9) + u_shape = 2 * (linear - 4) ** 2 + 17 + peak = -2 * (linear - 4) ** 2 + 55 + assert _count_turning_points(np.ones(9)) == 0 + assert _count_turning_points(linear) == 0 + assert _count_turning_points(u_shape) == 1 + assert _count_turning_points(peak) == 1 + + assert _count_turning_points(np.array([1, 2, 3, 4, 5, 4, 3, 2, 1])) == 1 + # Double peak is 3 points + assert _count_turning_points(np.array([1, 2, 3, 4, 2, 4, 3, 2, 1])) == 3 + # But only one if the dip isn't prominent + assert _count_turning_points(np.array([1, 2, 3, 4, 3.8, 4, 3, 2, 1])) == 1 diff --git a/webapp/src/components/tabContentComponents/backgroundDetectComponents/paneBackgroundDetect.vue b/webapp/src/components/tabContentComponents/backgroundDetectComponents/paneBackgroundDetect.vue index 4a41c055..9de2ec45 100644 --- a/webapp/src/components/tabContentComponents/backgroundDetectComponents/paneBackgroundDetect.vue +++ b/webapp/src/components/tabContentComponents/backgroundDetectComponents/paneBackgroundDetect.vue @@ -5,6 +5,9 @@
  • Configure
    +

    + {{ backgroundDetectorName }} +

    +