Calibrating a light path for retinal stimulation

This post describes using an image sensor to calibrate a projector used to stimulate retinal tissue in an electrophysiology experiment.

Designing a light path without testing it is relatively easy, as many of the setup's flaws will stay hidden. The experimenter's time and, more importantly, the life of the animal killed to perform an experiment can be wasted if a recording setup is not collecting the intended data. But as calibration is often tedious and will probably require custom hardware and software, it is not surprising when it is neglected. This post describes the hardware and software used to calibrate a light path that projects the image from a digital micromirror device (DMD) onto a flattened explanted retina.

We take measurements of the following:

1. homography between DMD and projected image
2. image sharpness
3. irradiance reaching the tissue, for each LED
4. an intensity map across the projected image for each LED

We also record the projection of a number of test images that are included in image sets used as checkpoints for the system.

To record images of the tissue plane, we position a naked image sensor in the place of the retinal tissue and take photos with the help of a Raspberry Pi. The GitHub project BadenLab/light-path-alignment contains the web server run on the Raspberry Pi.

Measuring alignment as a homography

We put a checkerboard image on the DMD and measure its projection on a camera sensor. These two images allow us to calculate the homography between the images sent to the projector and the actual image that lands on the target plane.

The homography estimated from the above pair, mapping stimulus pixels to camera pixels, is:

This transformation corresponds to:

• Scale. (scale x, scale y) = \((4.206, 4.210)\). One stimulus pixel covers 4.21 camera pixels, and the x and y magnifications agree to within 0.1%, so there is no stretch. The IMX519 has 1.22 µm square pixels, and the original image pixels correspond to 7.70 µm on the DMD, giving a physical magnification of (0.666, 0.667); in other words, the projected image is 66% the size of the DMD array in the projector, which has a (H, W) of (6.16 mm x 9.86 mm).
• Rotation: \(\theta = -0.49^\circ\). The projected image is rotated about half a degree relative to the sensor.
• Shear: \(k = 0.0002\), about 0.01°. Negligible.
• Translation: \(\matr{t} = (-369.4,\; 91.5)\). The stimulus is wider and shorter than the camera sensor. The stimulus origin lands 369 camera pixels to the left of the sensor's edge and 92 pixels down from its top, so the left of the stimulus is being cropped.
• Perspective: \(\matr{v} \approx (-6 \times 10^{-7},\; 9 \times 10^{-7})\). Effectively zero, meaning the DMD and sensor planes are parallel and there is no keystone.

Measuring this homography allows us to:

• align and center the projected image on the tissue
• know the magnification of the system and the size of a pixel on the tissue plane
• create the reverse mapping from the tissue plane to the stimulus space

The last point, creating the reverse mapping, is very powerful in that it allows us to correct for any misalignment and uneven intensity in software by manipulating our stimulus.

Below we use OpenCV to calculate a homography by passing in the source and target images. Although the problem reduces to a library call, there are some tedious details, such as creating the correct test image expected by the library (see OpenCV's guide to homography).

def estimate_homography(
    src_checkerboard_img, target_checkerboard_img, n_row_corners, n_col_corners
):
    total_corners = n_row_corners * n_col_corners
    src_gray = cv2.cvtColor(src_checkerboard_img, cv2.COLOR_RGB2GRAY)
    target_gray = cv2.cvtColor(target_checkerboard_img, cv2.COLOR_RGB2GRAY)
    # Normalize target to use the full [0, 255] range. Captures under
    # narrowband illumination (e.g. single LED) often have low contrast
    # which hurts corner detection.
    target_gray = cv2.normalize(target_gray, None, 0, 255, cv2.NORM_MINMAX)
    stim_corners = _detect_corners(src_gray, n_row_corners, n_col_corners)
    target_corners = _detect_corners(target_gray, n_row_corners, n_col_corners)

    pts_src = einops.rearrange(
        stim_corners, "n 1 c -> n c", n=total_corners, c=2
    )
    pts_dst = einops.rearrange(
        target_corners, "n 1 c -> n c", n=total_corners, c=2
    )
    pts_dst = _align_corner_ordering(pts_src, pts_dst)

    H, mask = cv2.findHomography(
        pts_src, pts_dst, method=cv2.RANSAC, ransacReprojThreshold=3.0
    )
    n_inliers = mask.sum()
    _logger.info(f"Homography inliers: {n_inliers}/{len(mask)}")

    # Decompose for diagnostics
    A = H[:2, :2] / H[2, 2]
    U, s, Vt = np.linalg.svd(A)
    R = U @ Vt
    if np.linalg.det(R) < 0:
        U[:, 1] *= -1
        R = U @ Vt
    rot_deg = np.degrees(np.arctan2(R[1, 0], R[0, 0]))
    _logger.debug(f"Scale factors: {s[0]:.4f}, {s[1]:.4f}")
    _logger.debug(f"Rotation: {rot_deg:.2f}°")
    _logger.debug(f"H:\n{H}")
    return H

When tweaking the XYZ stages in the light path, we run a live capture and stream it over HTTP to see quick feedback.

Image sharpness

The projected image can be well aligned while still being poorly focused. We estimate how well the light path is in focus by measuring the sharpness of the captured checkerboard image.

def sharpness(arr):
    """Laplacian variance focus measure.

    Args:
        arr: (H, W, C) uint8 array in [0, 255].

    Returns:
        Variance of the Laplacian (higher is sharper).
    """
    gray = cv2.cvtColor(arr, cv2.COLOR_RGB2GRAY)
    lap = cv2.Laplacian(gray, cv2.CV_64F)
    return float(lap.var())

Irradiance for each LED

We use a PM100D with a S130VC 200nm-1100nm optical power sensor placed where the tissue is to be mounted to measure the power reaching the tissue. We display a square on the DMD, small enough that it fits within the S130VC's sensor area. The power reading, paired with the magnification encoded in the measured homography and the known size of the square on the DMD allows us to calculate irradiance (power per area) for each LED connected to the projector. It is more valuable to record irradiance rather than simply power; if you record power only, the measurement is dependent on the magnification, square size and sensor area.

Taking the power meter readings is a bit time-consuming, so once they have been taken, we rely on diffs between test images to determine if LED intensities have drifted and need resetting.

Intensity map for each LED

Vignetting, illumination optics and device defects mean that the light path will not project a flat image with uniform intensity. We record an intensity map with the image sensor at the tissue plane, and we use this map to correct for uneven illumination.

Measuring an intensity map with dot grids

We present a grid of dots for each color at 16 different offsets, corresponding to 4 offset positions in each of the x and y directions. We combine the 16 captures to create an intensity map for each LED.

The three intensity maps:

If you are surprised to see the smudges, so were we! These smudges are a consequence of poor manufacturing by EKB, the company that makes the projectors. After ruling out defects with our sensor and any exposed surfaces, I have come to believe these smudges are on the surface of the prism that is placed against the DMD inside the projector. I placed a second EKB projector into the same light path and observed a different set of artifacts, including a clothing fiber.

Luckily, the retinal tissue is placed in the central area, taking up an 800x800 region of the 1280x800 display, which avoids most of the dirt artifacts.

In addition to the dirt artifacts, the lower frequency intensity pattern is different for each LED. We believe this is an issue with the design of EKB's illumination optics.

We used dot grids to measure the spatial intensity so that we were measuring impulse responses without the blur we expected if a full flat image was projected; however, later testing showed that flat images give similar intensity maps. As capturing a single image per LED is quicker and avoids having to merge multiple dot grid captures, I would recommend this approach first.

Intensity correction

By mapping the intensity measurements backwards through the estimated homography, we can place our intensity map in stimulus space. This allows us to modify an image to invert the effects of the light path, such as the uneven illumination.

Colorspace mapping

When we wish to present natural stimuli in the form of images or videos taken by cameras, we need an additional conversion to account for our display primaries.

For this, we need to run the conversion:

sRGB → linear RGB → XYZ → LED intensities → intensity-corrected LED intensities → gamma-encoded LED intensities

You can still see dark smudges caused by the unclean surfaces within the projector. Correcting for the uneven illumination eats up the usable intensity range, and for the above image, peak intensity was reduced by about 50% so that the darker areas could be driven harder to match. But this is still not enough to correct those smudges.