"""Functions to set up a Raspberry Pi Camera v2 for scientific use. This module provides slower, simpler functions to set the gain, exposure, and white balance of a Raspberry Pi camera, using the ``picamera2`` Python library. It's mostly used by the OpenFlexure Microscope, though it deliberately has no hard dependencies on said software, so that it's useful on its own. There are three main calibration steps: * Setting exposure time and gain to get a reasonably bright image. * Fixing the white balance to get a neutral image * Taking a uniform white image and using it to calibrate the Lens Shading Table The most reliable way to do this, avoiding any issues relating to "memory" or nonlinearities in the camera's image processing pipeline, is to use raw images. This is quite slow, but very reliable. The three steps above can be accomplished by: .. code-block:: python picamera = picamera2.Picamera2() adjust_shutter_and_gain_from_raw(picamera) adjust_white_balance_from_raw(picamera) lst = lst_from_camera(picamera) picamera.lens_shading_table = lst """ # Disable N806 & 803, which checks that all variables and args are lowercase. # This is due to the number of matrix calculations and colour channel # calculations that are clearer using the standard R, G, B, or L, Cr, Cb terms. # ruff: noqa: N806 N803 from __future__ import annotations import gc import logging import time from typing import List, Literal, Optional, Tuple from pydantic import BaseModel import numpy as np from scipy.ndimage import zoom from picamera2 import Picamera2 import picamera2 LensShadingTables = tuple[np.ndarray, np.ndarray, np.ndarray] def load_default_tuning(cam: Picamera2) -> dict: """Load the default tuning file for the camera. This will open and close the camera to determine its model. If you are using a model that's supported by ``picamera2`` it should have a tuning file built in. If not, this will probably crash with an error. Error handling for unsupported cameras is not something we are likely to test in the short term. """ cp = cam.camera_properties fname = f"{cp['Model']}.json" try: return cam.load_tuning_file(fname) except RuntimeError: tuning_dir = "/usr/share/libcamera/ipa/raspberrypi" # from picamera2 v0.3.9 # The directory above has been removed from the search path seems # odd - as that's where the files currently are on a default # Raspbian image. This may need updating if the files have moved # in future updates to the system libcamera package return cam.load_tuning_file(fname, dir=tuning_dir) def set_minimum_exposure(camera: Picamera2) -> None: """Enable manual exposure, with low gain and shutter speed. We set exposure mode to manual, analog and digital gain to 1, and shutter speed to the minimum (8us for Pi Camera v2) Note ISO is left at auto, because this is needed for the gains to be set correctly. """ # Disable Automatic exposure and gain algorithm (AeEnable), and set analogue # gain and exposure time. # Setting the shutter speed to 1us will result in it being set # to the minimum possible, which is ~8us for PiCamera v2 camera.set_controls({"AeEnable": False, "AnalogueGain": 1, "ExposureTime": 1}) time.sleep(0.5) class ExposureTest(BaseModel): """Record the results of testing the camera's current exposure settings.""" level: int exposure_time: int analog_gain: float def test_exposure_settings(camera: Picamera2, percentile: float) -> ExposureTest: """Evaluate current exposure settings using a raw image. CAMERA SHOULD BE STARTED! We will acquire a raw image and calculate the given percentile of the pixel values. We return a dictionary containing the percentile (which will be compared to the target), as well as the camera's shutter and gain values. """ camera.capture_array("raw") # controls might not be updated for the first frame? max_brightness = np.percentile( channels_from_bayer_array(camera.capture_array("raw")), percentile, ) # The reported brightness can, theoretically, be negative or zero # because of black level compensation. The line below forces a # minimum value of 1 which will keep things well-behaved! if max_brightness < 1: logging.warning( f"Measured brightness of {max_brightness}. " "This should normally be >= 1, and may indicate the " "camera's black level compensation has gone wrong." ) max_brightness = 1 metadata = camera.capture_metadata() result = ExposureTest( level=max_brightness, exposure_time=int(metadata["ExposureTime"]), analog_gain=float(metadata["AnalogueGain"]), ) logging.info(f"{result.model_dump()}") return result def check_convergence(test: ExposureTest, target: int, tolerance: float) -> bool: """Check whether the brightness is within the specified target range.""" return abs(test.level - target) < target * tolerance def adjust_shutter_and_gain_from_raw( camera: Picamera2, target_white_level: int = 700, max_iterations: int = 20, tolerance: float = 0.05, percentile: float = 99.9, ) -> float: """Adjust exposure and analog gain based on raw images. This routine is slow but effective. It uses raw images, so we are not affected by white balance or digital gain. :param camera: A Picamera2 object. :param target_white_level: The raw, 10-bit value we aim for. The brightest pixels should be approximately this bright. Maximum possible is about 900, 700 is reasonable. :param max_iterations: We will terminate once we perform this many iterations, whether or not we converge. More than 10 shouldn't happen. :param tolerance: How close to the target value we consider "done". Expressed as a fraction of the ``target_white_level`` so 0.05 means +/- 5% :param percentile: Rather then use the maximum value for each channel, we calculate a percentile. This makes us robust to single pixels that are bright/noisy. 99.9% still picks the top of the brightness range, but seems much more reliable than just ``np.max()``. """ # TODO: read black level and bit depth from camera? if target_white_level * (tolerance + 1) >= 959: raise ValueError( "The target level is too high - a saturated image would be " "considered successful. target_white_level * (tolerance + 1) " "must be less than 959." ) config = camera.create_still_configuration(raw={"format": "SBGGR10"}) camera.configure(config) camera.start() set_minimum_exposure(camera) # We start with very low exposure settings and work up # until either the brightness is high enough, or we can't increase the # shutter speed any more. iterations = 0 while iterations < max_iterations: test = test_exposure_settings(camera, percentile) if check_convergence(test, target_white_level, tolerance): break iterations += 1 # Adjust shutter speed so that the brightness approximates the target # NB we put a maximum of 8 on this, to stop it increasing too quickly. new_time = int(test.exposure_time * min(target_white_level / test.level, 8)) camera.controls.ExposureTime = new_time camera.controls.AeEnable = False time.sleep(0.5) # Check whether the shutter speed is still going up - if not, we've hit a maximum if camera.capture_metadata()["ExposureTime"] == test.exposure_time: logging.info(f"Shutter speed has maxed out at {test.exposure_time}") break # Now, if we've not converged, increase gain until we converge or run out of options. while iterations < max_iterations: test = test_exposure_settings(camera, percentile) if check_convergence(test, target_white_level, tolerance): break iterations += 1 # Adjust gain to make the white level hit the target, again with a maximum camera.controls.AnalogueGain = test.analog_gain * min( target_white_level / test.level, 2 ) time.sleep(0.5) # Check the gain is still changing - if not, we have probably hit the maximum if camera.capture_metadata()["AnalogueGain"] == test.analog_gain: logging.info(f"Gain has maxed out at {test.analog_gain}") break if check_convergence(test, target_white_level, tolerance): logging.info(f"Brightness has converged to within {tolerance * 100:.0f}%.") else: logging.warning( f"Failed to reach target brightness of {target_white_level}." f"Brightness reached {test.level} after {iterations} iterations." ) return test.level def adjust_white_balance_from_raw( camera: Picamera2, percentile: float = 99, luminance: Optional[np.ndarray] = None, Cr: Optional[np.ndarray] = None, Cb: Optional[np.ndarray] = None, luminance_power: float = 1.0, method: Literal["percentile", "centre"] = "centre", ) -> Tuple[float, float]: """Adjust the white balance in a single shot, based on the raw image. NB if ``channels_from_raw_image`` is broken, this will go haywire. We should probably have better logic to verify the channels really are BGGR... """ config = camera.create_still_configuration(raw={"format": "SBGGR10"}) camera.configure(config) camera.start() channels = channels_from_bayer_array(camera.capture_array("raw")) # TODO: read black level from camera rather than hard-coding 64 blacklevel = 64 if luminance is not None and Cr is not None and Cb is not None: # Reconstruct a low-resolution image from the lens shading tables # and use it to normalise the raw image, to compensate for # the brightest pixels in each channel not coinciding. grids = grids_from_lst(np.array(luminance) ** luminance_power, Cr, Cb) channel_gains = 1 / grids if channel_gains.shape[1:] != channels.shape[1:]: channel_gains = upsample_channels(channel_gains, channels.shape[1:]) logging.info( f"Before gains, channel maxima are {np.max(channels, axis=(1, 2))}" ) channels = channels * channel_gains logging.info(f"After gains, channel maxima are {np.max(channels, axis=(1, 2))}") if method == "centre": _, height, width = channels.shape # Cut out the central 10% from 9/20 to 11/20... low_y_range = 9 * height // 20 hi_y_range = 11 * height // 20 low_x_range = 9 * width // 20 hi_x_range = 11 * width // 20 # ... and then take the mean of each bayer channel. centre_means = np.mean( channels[:, low_y_range:hi_y_range, low_x_range:hi_x_range], axis=(1, 2), ) # Subtract blacklevel before splitting into channels blue, g1, g2, red = centre_means - blacklevel else: blue, g1, g2, red = ( np.percentile(channels, percentile, axis=(1, 2)) - blacklevel ) green = (g1 + g2) / 2.0 new_awb_gains = (green / red, green / blue) if Cr is not None and Cb is not None: # The LST algorithm normalises Cr and Cb by their minimum. # The lens shading correction only ever boosts the red and blue values. # Here, we decrease the gains by the minimum value of Cr and Cb. new_awb_gains = (green / red * np.min(Cr), green / blue * np.min(Cb)) logging.info( f"Raw white point is R: {red} G: {green} B: {blue}, " f"setting AWB gains to ({new_awb_gains[0]:.2f}, " f"{new_awb_gains[1]:.2f})." ) camera.controls.AwbEnable = False camera.controls.ColourGains = new_awb_gains time.sleep(0.2) m = camera.capture_metadata() print(f"Camera confirms gains are now {m['ColourGains']}") return new_awb_gains def channels_from_bayer_array(bayer_array: np.ndarray) -> np.ndarray: """Given the 'array' from a PiBayerArray, return the 4 channels.""" bayer_pattern: List[Tuple[int, int]] = [(0, 0), (0, 1), (1, 0), (1, 1)] bayer_array = bayer_array.view(np.uint16) channels_shape: Tuple[int, int, int] = ( 4, bayer_array.shape[0] // 2, bayer_array.shape[1] // 2, ) channels: np.ndarray = np.zeros(channels_shape, dtype=bayer_array.dtype) for i, offset in enumerate(bayer_pattern): # We simplify life by dealing with only one channel at a time. channels[i, :, :] = bayer_array[offset[0] :: 2, offset[1] :: 2] return channels def get_16x12_grid(chan: np.ndarray, dx: int, dy: int) -> np.ndarray: """Compresses channel down to a 16x12 grid - from libcamera. This is taken from https://git.linuxtv.org/libcamera.git/tree/utils/raspberrypi/ctt/ctt_alsc.py for consistency. """ grid = [] # since left and bottom border will not necessarily have rectangles of # dimension dx x dy, the final iteration has to be handled separately. for i in range(11): for j in range(15): grid.append(np.mean(chan[dy * i : dy * (1 + i), dx * j : dx * (1 + j)])) grid.append(np.mean(chan[dy * i : dy * (1 + i), 15 * dx :])) for j in range(15): grid.append(np.mean(chan[11 * dy :, dx * j : dx * (1 + j)])) grid.append(np.mean(chan[11 * dy :, 15 * dx :])) # return as np.array, ready for further manipulation return np.reshape(np.array(grid), (12, 16)) def upsample_channels(grids: np.ndarray, shape: tuple[int]) -> np.ndarray: """Zoom an image in the last two dimensions. This is effectively the inverse operation of ``get_16x12_grid`` """ zoom_factors = [ 1, ] + list(np.ceil(np.array(shape) / np.array(grids.shape[1:]))) return zoom(grids, zoom_factors, order=1)[:, : shape[0], : shape[1]] def downsampled_channels(channels: np.ndarray, blacklevel=64) -> list[np.ndarray]: """Generate a downsampled, un-normalised image from which to calculate the LST. TODO: blacklevel probably ought to be determined from the camera... """ channel_shape = np.array(channels.shape[1:]) lst_shape = np.array([12, 16]) step = np.ceil(channel_shape / lst_shape).astype(int) return np.stack( [ get_16x12_grid( channels[i, ...].astype(float) - blacklevel, step[1], step[0] ) for i in range(channels.shape[0]) ], axis=0, ) def lst_from_channels(channels: np.ndarray) -> LensShadingTables: """Given the 4 Bayer colour channels from a white image, generate a LST. Internally, is just calls ``downsampled_channels`` and ``lst_from_grids``. """ grids = downsampled_channels(channels) return lst_from_grids(grids) def lst_from_grids(grids: np.ndarray) -> LensShadingTables: """Given 4 downsampled grids, generate the luminance and chrominance tables. The grids are the 4 BAYER channels RGGB The LST format has changed with ``picamera2`` and now uses a fixed resolution, and is in luminance, Cr, Cb format. This function returns three ndarrays of luminance, Cr, Cb, each with shape (12, 16). """ # Calculated red, green, and blue channels from Bayer data r: np.ndarray = grids[3, ...] g: np.ndarray = np.mean(grids[1:3, ...], axis=0) b: np.ndarray = grids[0, ...] # What we actually want to calculate is the gains needed to compensate for the # lens shading - that's 1/lens_shading_table_float as we currently have it. # Minimum luminance gain is 1 luminance_gains: np.ndarray = np.max(g) / g cr_gains: np.ndarray = g / r cb_gains: np.ndarray = g / b return luminance_gains, cr_gains, cb_gains def grids_from_lst(lum: np.ndarray, Cr: np.ndarray, Cb: np.ndarray) -> np.ndarray: """Convert form luminance/chrominance dict to four RGGB channels. Note that these will be normalised - the maximum green value is always 1. Also, note that the channels are BGGR, to be consistent with the ``channels_from_raw_image`` function. This should probably change in the future. """ G = 1 / np.array(lum) R = G / np.array(Cr) B = G / np.array(Cb) return np.stack([B, G, G, R], axis=0) def set_static_lst( tuning: dict, luminance: np.ndarray, cr: np.ndarray, cb: np.ndarray, ) -> None: """Update the ``rpi.alsc`` section of a camera tuning dict to use a static correction. ``tuning`` will be updated in-place to set its shading to static, and disable any adaptive tweaking by the algorithm. """ for table in luminance, cr, cb: assert np.array(table).shape == (12, 16), "Lens shading tables must be 12x16!" alsc = Picamera2.find_tuning_algo(tuning, "rpi.alsc") alsc["n_iter"] = 0 # disable the adaptive part alsc["luminance_strength"] = 1.0 alsc["calibrations_Cr"] = [ {"ct": 4500, "table": as_flat_rounded_list(cr, round_to=3)} ] alsc["calibrations_Cb"] = [ {"ct": 4500, "table": as_flat_rounded_list(cb, round_to=3)} ] alsc["luminance_lut"] = as_flat_rounded_list(luminance, round_to=3) def set_static_ccm( tuning: dict, col_corr_matrix: tuple[ float, float, float, float, float, float, float, float, float ], ) -> None: """Update the ``rpi.alsc`` section of a camera tuning dict to use a static correction. ``tuning`` will be updated in-place to set its shading to static, and disable any adaptive tweaking by the algorithm. """ ccm = Picamera2.find_tuning_algo(tuning, "rpi.ccm") ccm["ccms"] = [{"ct": 2860, "ccm": col_corr_matrix}] def get_static_ccm(tuning: dict) -> None: """Get the ``rpi.ccm`` section of a camera tuning dict.""" ccm = Picamera2.find_tuning_algo(tuning, "rpi.ccm") return ccm["ccms"] def lst_is_static(tuning: dict) -> bool: """Whether the lens shading table is set to static.""" alsc = Picamera2.find_tuning_algo(tuning, "rpi.alsc") return alsc["n_iter"] == 0 def set_static_geq( tuning: dict, offset: int = 65535, ) -> None: """Update the ``rpi.geq`` section of a camera tuning dict. :param tuning: the raspberry pi tuning file. This will be updated in-place to set the geq offset to the given value. :param offset: The desired green equalisation offset. Default 65535. The default is the maximum allowed value. This means the brightness will always be below the threshold where averaging is used. This is default as we always need the green equalisation to averages the green pixels in the red and blue rows due to the chief ray angle compensation issue when the the stock lens is replaced by an objective. """ geq = Picamera2.find_tuning_algo(tuning, "rpi.geq") geq["offset"] = offset # max out offset to disable the adaptive green equalisation def _geq_is_static(tuning: dict) -> bool: """Whether the green equalisation is set to static.""" geq = Picamera2.find_tuning_algo(tuning, "rpi.geq") return geq["offset"] == 65535 def index_of_algorithm(algorithms: list[dict], algorithm: str) -> int: """Find the index of an algorithm's section in the tuning file.""" for i, a in enumerate(algorithms): if algorithm in a: return i raise ValueError(f"Algorithm {algorithm} is not available.") def copy_alsc_section(from_tuning: dict, to_tuning: dict) -> None: """Copy the ``rpi.alsc`` algorithm from one tuning to another. This is done in-place, i.e. modifying to_tuning. """ # Using Picamera2 function to find the relevant sub-dict for each tuning file from_i = index_of_algorithm(from_tuning["algorithms"], "rpi.alsc") to_i = index_of_algorithm(to_tuning["algorithms"], "rpi.alsc") # Updating the dictionary in place. to_tuning["algorithms"][to_i] = from_tuning["algorithms"][from_i] def lst_from_camera(camera: Picamera2) -> LensShadingTables: """Acquire a raw image and use it to calculate a lens shading table.""" channels = raw_channels_from_camera(camera) return lst_from_channels(channels) def raw_channels_from_camera(camera: Picamera2) -> LensShadingTables: """Acquire a raw image and return a 4xNxM array of the colour channels.""" if camera.started: camera.stop_recording() # We will acquire a raw image with unpacked pixels, which is what the # format below requests. Bit depth and Bayer order may be overwritten. # TODO: don't assume 10-bit - the high quality camera uses 12. # TODO: what's the best mode to use here? config = camera.create_still_configuration(raw={"format": "SBGGR10"}) camera.configure(config) camera.start() raw_image = camera.capture_array("raw") camera.stop() # Now we need to calculate a lens shading table that would make this flat. # raw_image is a 3D array, with full resolution and 3 colour channels. No # de-mosaicing has been done, so 2/3 of the values are zero (3/4 for R and B # channels, 1/2 for green because there's twice as many green pixels). raw_format = camera.camera_configuration()["raw"]["format"] print(f"Acquired a raw image in format {raw_format}") return channels_from_bayer_array(raw_image) def recreate_camera_manager() -> None: """Delete and recreate the camera manager. This is necessary to ensure the tuning file is re-read. """ del Picamera2._cm gc.collect() Picamera2._cm = picamera2.picamera2.CameraManager() def as_flat_rounded_list(array: np.ndarray, round_to: int = 3) -> list[float]: """Flatten array, round, and then convert to list.""" return np.reshape(array, -1).round(round_to).tolist()