Skip to content

Convolve with a filter bank in peak finding #70

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
HazenBabcock merged 3 commits into from
Sep 13, 2026
Merged

Convolve with a filter bank in peak finding #70

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
564dee3
Add a filter bank to matched_filter.c
HazenBabcock Sep 12, 2026
8b98701
Use a filter bank in the spliner peak finder
HazenBabcock Sep 12, 2026
eec2450
Use filter banks in the multiplane peak finder
HazenBabcock Sep 12, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
112 changes: 62 additions & 50 deletions storm_analysis/multi_plane/fitting_mp.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,10 +57,11 @@ def __init__(self, parameters = None, **kwds):
self.backgrounds = []
self.bg_filters = [] # Matched filters for background estimation.
self.images = []
self.mfilters = [] # Matched filters for foreground estimation.
self.mfilter_banks = [] # Foreground matched filter banks, one per channel.
self.n_channels = None
self.n_zvals = 0 # The number of z values, i.e. the size of each bank.
self.variances = [] # sCMOS variance data.
self.vfilters = [] # Matched filters for (background) variance estimation.
self.vfilter_banks = [] # (Background) variance matched filter banks, one per channel.
self.xt = []
self.yt = []
self.z_values = []
Expand All @@ -77,10 +78,9 @@ def cleanUp(self):
at.cleanup()

# Clean up foreground filters and variance filters.
for i in range(len(self.mfilters)):
for j in range(len(self.mfilters[i])):
self.mfilters[i][j].cleanup()
self.vfilters[i][j].cleanup()
for i in range(len(self.mfilter_banks)):
self.mfilter_banks[i].cleanup()
self.vfilter_banks[i].cleanup()

# Clean up background filters.
for i in range(len(self.bg_filters)):
Expand Down Expand Up @@ -197,19 +197,26 @@ def peakFinder(self, fit_peaks_images):
for fi in fit_peaks_images:
tf.write(fi.astype(numpy.float32))

#
# Each channel's image is the same at every z value, so convolve it
# with all of that channel's variance filters in one pass rather than
# transforming it again for each z value.
#
# I believe that this is correct, the variance of the weighted average
# of independent processes is calculated using the square of the weights.
#
conv_vars = []
for j in range(self.n_channels):
conv_vars.append(self.vfilter_banks[j].convolve(fit_peaks_images[j] + self.backgrounds[j]))

# Iterate over z values.
for i in range(len(self.vfilters)):
for i in range(self.n_zvals):
bg_variance = numpy.zeros(fit_peaks_images[0].shape)

# Iterate over channels / planes.
for j in range(len(self.vfilters[i])):
for j in range(self.n_channels):

# Convolve fit image + background with the appropriate variance filter.
#
# I believe that this is correct, the variance of the weighted average
# of independent processes is calculated using the square of the weights.
#
conv_var = self.vfilters[i][j].convolve(fit_peaks_images[j] + self.backgrounds[j])
conv_var = conv_vars[j][i]

# Transform variance to the channel 0 frame.
if self.atrans[j] is None:
Expand Down Expand Up @@ -240,16 +247,20 @@ def peakFinder(self, fit_peaks_images):
fg_averages = [] # This is the average foreground across all the planes for each z value.
foregrounds = [] # This is the foreground for each plane and z value.

# As above, one pass per channel rather than one per (z, channel).
conv_fgs = []
for j in range(self.n_channels):
conv_fgs.append(self.mfilter_banks[j].convolve(self.images[j] - fit_peaks_images[j] - self.backgrounds[j]))

# Iterate over z values.
for i in range(len(self.mfilters)):
for i in range(self.n_zvals):
foreground = numpy.zeros(fit_peaks_images[0].shape)
foregrounds.append([])

# Iterate over channels / planes.
for j in range(len(self.mfilters[i])):
for j in range(self.n_channels):

# Convolve image / background with the appropriate PSF.
conv_fg = self.mfilters[i][j].convolve(self.images[j] - fit_peaks_images[j] - self.backgrounds[j])
conv_fg = conv_fgs[j][i]

# Store convolved image in foregrounds.
foregrounds[i].append(conv_fg)
Expand Down Expand Up @@ -295,7 +306,7 @@ def setVariances(self, variances):
We initialize the following here because at __init__ we
don't know how big the images are. This is called after
the analysis specific version pads the variances and
creates the mfilters[] and the vfilters[] class members.
creates the mfilter_banks[] and the vfilter_banks[] class members.

Note the assumption that every frame in all the movies
is the same size.
Expand Down Expand Up @@ -332,16 +343,21 @@ def setVariances(self, variances):
# we are applying to the foreground.
#

# One pass per channel, as in peakFinder().
#
conv_vars = []
for j in range(self.n_channels):
conv_vars.append(self.vfilter_banks[j].convolve(variances[j]))

# Iterate over z values.
#
for i in range(len(self.mfilters)):
for i in range(self.n_zvals):
variance = numpy.zeros(variances[0].shape)

# Iterate over channels / planes.
for j in range(len(self.mfilters[i])):
for j in range(self.n_channels):

# Convolve variance with the appropriate variance filter.
conv_var = self.vfilters[i][j].convolve(variances[j])
conv_var = conv_vars[j][i]

# Transform variance to the channel 0 frame.
if self.atrans[j] is None:
Expand Down Expand Up @@ -451,13 +467,14 @@ def setVariances(self, variances):

# Create "foreground" and "variance" filters.
#
# These are stored in a list indexed by z value, then by
# channel / plane. So self.mfilters[1][2] is the filter
# for z value 1, plane 2.
# These are stored as one bank per channel / plane, each bank
# holding that channel's filter for every z value. So
# self.mfilter_banks[2] convolves plane 2 at all z values at once.
#
mfilter_psfs = [[] for i in range(self.n_channels)]
vfilter_psfs = [[] for i in range(self.n_channels)]

for i, mfilter_z in enumerate(self.mfilters_z):
self.mfilters.append([])
self.vfilters.append([])

for j, psf_object in enumerate(self.psf_objects):
psf = psf_object.getPSF(mfilter_z,
Expand All @@ -469,21 +486,22 @@ def setVariances(self, variances):
# or if it does that they are very small.
#
psf_norm = psf/numpy.sum(psf)
self.mfilters[i].append(matchedFilterC.MatchedFilter(psf_norm,
fftw_estimate = self.parameters.getAttr("fftw_estimate"),
memoize = True,
max_diff = 1.0e-3))
self.vfilters[i].append(matchedFilterC.MatchedFilter(psf_norm * psf_norm,
fftw_estimate = self.parameters.getAttr("fftw_estimate"),
memoize = True,
max_diff = 1.0e-3))
mfilter_psfs[j].append(psf_norm)
vfilter_psfs[j].append(psf_norm * psf_norm)

# Save a pictures of the PSFs for debugging purposes.
if self.check_mode:
print("psf max", numpy.max(psf))
filename = "psf_z{0:.3f}_c{1:d}.tif".format(mfilter_z, j)
tifffile.imwrite(filename, psf.astype(numpy.float32))

self.n_zvals = len(self.mfilters_z)
for j in range(self.n_channels):
self.mfilter_banks.append(matchedFilterC.MatchedFilterBank(mfilter_psfs[j],
fftw_estimate = self.parameters.getAttr("fftw_estimate")))
self.vfilter_banks.append(matchedFilterC.MatchedFilterBank(vfilter_psfs[j],
fftw_estimate = self.parameters.getAttr("fftw_estimate")))

# This handles the rest of the initialization.
#
super(MPPeakFinderArb, self).setVariances(variances)
Expand Down Expand Up @@ -578,25 +596,19 @@ def setVariances(self, variances):
# Create "foreground" and "variance" filters. There is
# only one z value here.
#
# These are stored in a list indexed by z value, then by
# channel / plane. So self.mfilters[1][2] is the filter
# for z value 1, plane 2.
# These are stored as one bank per channel / plane, each bank
# holding that channel's filter for every z value. So
# self.mfilter_banks[2] convolves plane 2 at all z values at once.
#
self.mfilters.append([])
self.vfilters.append([])

psf_norm = fitting.gaussianPSF(variances[0].shape, self.parameters.getAttr("foreground_sigma"))
var_norm = psf_norm * psf_norm

self.n_zvals = 1
for i in range(self.n_channels):
self.mfilters[0].append(matchedFilterC.MatchedFilter(psf_norm,
fftw_estimate = self.parameters.getAttr("fftw_estimate"),
memoize = True,
max_diff = 1.0e-3))
self.vfilters[0].append(matchedFilterC.MatchedFilter(var_norm,
fftw_estimate = self.parameters.getAttr("fftw_estimate"),
memoize = True,
max_diff = 1.0e-3))
self.mfilter_banks.append(matchedFilterC.MatchedFilterBank([psf_norm],
fftw_estimate = self.parameters.getAttr("fftw_estimate")))
self.vfilter_banks.append(matchedFilterC.MatchedFilterBank([var_norm],
fftw_estimate = self.parameters.getAttr("fftw_estimate")))

# Save a pictures of the PSFs for debugging purposes.
if self.check_mode:
Expand Down
49 changes: 30 additions & 19 deletions storm_analysis/sa_library/fitting.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -619,9 +619,9 @@ def __init__(self, parameters = None, psf_object = None, **kwds):
kwds["parameters"] = parameters
super(PeakFinderArbitraryPSF, self).__init__(**kwds)

self.fg_mfilter = []
self.fg_mfilter = None # Foreground MatchedFilterBank, one psf per z value.
self.fg_mfilter_zval = []
self.fg_vfilter = []
self.fg_vfilter = None # Foreground variance MatchedFilterBank.
self.psf_object = psf_object
self.z_values = []

Expand Down Expand Up @@ -666,9 +666,9 @@ def __init__(self, parameters = None, psf_object = None, **kwds):

def cleanUp(self):
super(PeakFinderArbitraryPSF, self).cleanUp()
for i in range(len(self.fg_mfilter)):
self.fg_mfilter[i].cleanup()
self.fg_vfilter[i].cleanup()
if self.fg_mfilter is not None:
self.fg_mfilter.cleanup()
self.fg_vfilter.cleanup()

def newImage(self, new_image):
"""
Expand All @@ -682,27 +682,32 @@ def newImage(self, new_image):
# If does not already exist, create filter objects from
# the PSF at different z values.
#
if (len(self.fg_mfilter) == 0):
if self.fg_mfilter is None:
mfilter_psfs = []
vfilter_psfs = []
for zval in self.fg_mfilter_zval:
psf = self.psf_object.getPSF(zval,
shape = new_image.shape,
normalize = False)
psf_norm = psf/numpy.sum(psf)
fg_mfilter = matchedFilterC.MatchedFilter(psf_norm,
fftw_estimate = self.parameters.getAttr("fftw_estimate"),
memoize = True,
max_diff = 1.0e-3)
self.fg_mfilter.append(fg_mfilter)
self.fg_vfilter.append(matchedFilterC.MatchedFilter(psf_norm * psf_norm,
fftw_estimate = self.parameters.getAttr("fftw_estimate"),
memoize = True,
max_diff = 1.0e-3))
mfilter_psfs.append(psf_norm)
vfilter_psfs.append(psf_norm * psf_norm)

# Save a picture of the PSF for debugging purposes.
if self.check_mode:
print("psf max", numpy.max(psf))
filename = "psf_{0:.3f}.tif".format(zval)
tifffile.imwrite(filename, psf.astype(numpy.float32))

#
# These are banks rather than one filter per z value because
# peakFinder() applies all of them to the same image, so the
# forward FFT only needs to happen once per bank.
#
self.fg_mfilter = matchedFilterC.MatchedFilterBank(mfilter_psfs,
fftw_estimate = self.parameters.getAttr("fftw_estimate"))
self.fg_vfilter = matchedFilterC.MatchedFilterBank(vfilter_psfs,
fftw_estimate = self.parameters.getAttr("fftw_estimate"))

def peakFinder(self, fit_peaks_image):
"""
Expand Down Expand Up @@ -736,11 +741,18 @@ def peakFinder(self, fit_peaks_image):
fg_tif = tifffile.TiffWriter("foreground.tif")
fg_bg_ratio_tif = tifffile.TiffWriter("fg_bg_ratio.tif")

#
# Both of these inputs are the same at every z value, so each bank
# transforms its input once instead of once per z value.
#
backgrounds = self.fg_vfilter.convolve(bg_var)
foregrounds = self.fg_mfilter.convolve(self.image - self.background - fit_peaks_image)

masked_images = []
for i in range(len(self.fg_mfilter)):
for i in range(len(self.fg_mfilter_zval)):

# Estimate background variance at this particular z value.
background = self.fg_vfilter[i].convolve(bg_var)
background = backgrounds[i]

# Remove problematic values.
#
Expand All @@ -754,8 +766,7 @@ def peakFinder(self, fit_peaks_image):
bg_std = numpy.sqrt(background)

# Calculate foreground.
foreground = self.image - self.background - fit_peaks_image
foreground = self.fg_mfilter[i].convolve(foreground)
foreground = foregrounds[i]

# Calculate foreground in units of signal to noise.
fg_bg_ratio = foreground/bg_std
Expand Down
Loading