Find focus by optimizing a transverse spatial frequency band

[talonchandler]

This PR implement a focus_from_transverse_band function that (by default) finds a focal plane by maximizing the transverse mid-band power. By inspecting the power spectrum of brightfield data from cells (from both Altos and hummingbird), I defined the midband as betwen 1/8 and 1/4 of the incoherent cutoff frequency (2*NA/lambda).

Usage: Basic usage with default options looks like:

>>> zyx_array.shape
(11, 2048, 2048)
>>> from waveorder.focus import focus_from_transverse_band
>>> slice = focus_from_transverse_band(zyx_array, NA_det=0.55, lambda_ill=0.532, pixel_size=6.5/20)
>>> in_focus_data = data[slice,:,:]

The function also has two optional arguments:

• midband_fractions allows the user to use an arbitrary definition of the "midband"
• mode allows the user to choose between maximizing or minimizing the power in this band.

I've chosen the default arguments empirically based on label-free brightfield data, but future users may need to find the focus from fluorescence or other modalities and these options are intended to be a good starting place.

Please see the docstring for more details.

Testing on a dataset: Here's a complete example snippet to get started:

import napari
v = napari.Viewer()
from waveorder.io import WaveorderReader
from waveorder.focus import focus_from_transverese_band

reader = WaveorderReader(`/hpc/projects/compmicro/rawdata/hummingbird/Janie/2022_03_15_orgs_nuc_mem_63x_04NA/all_21_3.zarr')
bf_idx = reader.channel_names.index('BF')
for i in range(reader.get_num_positions()):
    arr = reader.get_array(i)[0,bf_idx,...]
    focus_idx = focus_from_transverse_band(arr, NA_det=0.55, lambda_ill=0.532, pixel_size=6.5/20)
    print(f'Position {i}, focus slice = {focus_idx}')
    v.add_image(arr)
    v.dims.set_point(0, focus_idx)

For this dataset, the in-focus slices for the first few positions are: 26, 70, 64, 66, 55, 54.

Limitations

• Model-free: this algorithm does not use an imaging model...future iterations might be integrated with knowledge of the full PSFs/OTFs for better performance.
• Speed: this algorithm takes a 2D FFT for each slice (more expensive than a 3D FFT), which can be time consuming: a 97 x 2048 x 2048 volumes (roughly our largest single volume) takes about 40 seconds. If speed becomes a bottleneck there are many improvements that can be made.

Next steps: @Soorya19Pradeep I understand that you're interested in applying this immediately. Can you try it on your datasets and report back? Send me a slack if you're unfamiliar with running code from an unmerged branch like this.

@ziw-liu I'd appreciate your code review.

[talonchandler]

ImportError on typing in 3.7. @ziw-liu do you think we're ready to remove 3.7 support?

[ziw-liu]

Sorry that I didn't mention that PEP-586 needs 3.8. We can still work around this if 3.7 is absolutely needed.

do you think we're ready to remove 3.7 support?

Theoretically yes. NumPy has dropped support for it and the Python version itself will go EOL in months. We are fine to move on (and start picking up 3.11) if there is no user workflow that has to be in 3.7.

Edit: just noting that a separate issue/PR may be appropriate if we were to make the move

[Soorya19Pradeep]

I have tested the script for finding focus in phase image stacks and fluorescence channel stacks (nucleus and membrane) in Vero cells at low and high confluence. I have also tested the script on very low signal-to-noise ratio dataset, both phase images and membrane channels. All these datasets were handled perfectly for focus finding by the script. I have also tested the script on raw brightfield images of Vero cells, which worked as well.