From 564dee31c380ef46455a3b8aa9cf6a7ef4c4b7a6 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Sat, 12 Sep 2026 19:09:39 -0400 Subject: [PATCH 1/3] Add a filter bank to matched_filter.c Applying N psfs to one image with N filter structures computes the forward FFT of that image N times, and all N results are identical. A bank keeps one pair of FFTW plans and N precomputed psf transforms, transforms the image once, and then does N multiply-and-invert passes. The transform count per call goes from 2N to N+1. image_fft is separate from work_fft because convolve() does its complex multiply in place. Reusing one buffer would destroy the very transform the bank exists to share. Measured on 270x270, the size the 256x256 test data pads to, interleaving the two arms so that neither gets a whole block of warm cache: N separate bank speedup 2N/(N+1) 1 0.568 ms 0.574 ms 0.99x 1.00x 2 1.267 0.935 1.36x 1.33x 3 1.911 1.317 1.45x 1.50x 5 3.339 1.997 1.67x 1.67x 8 5.499 3.106 1.77x 1.78x So it tracks the predicted ratio, and at N = 1 it correctly does nothing. There is no point using a bank for a single psf. The bank does not memoize. convolveMemo() only helps when the same filter sees nearly the same image twice in a row, and instrumenting the spliner peak finder showed that never happening once in 220 calls - the image changes every iteration as peaks are subtracted. Nothing is given up by leaving it out. test_matched_filter_bank1 asserts array_equal against N separate filters rather than allclose, since the bank changes only where the forward transform happens, not the arithmetic. Also lists the existing tests in the __main__ block, which called only test_matched_filter5. --- storm_analysis/sa_library/matched_filter.c | 152 ++++++++++++++++++ storm_analysis/sa_library/matched_filter_c.py | 68 ++++++++ storm_analysis/test/test_matched_filter.py | 82 ++++++++++ 3 files changed, 302 insertions(+) diff --git a/storm_analysis/sa_library/matched_filter.c b/storm_analysis/sa_library/matched_filter.c index 8469c0a8..742c9226 100644 --- a/storm_analysis/sa_library/matched_filter.c +++ b/storm_analysis/sa_library/matched_filter.c @@ -49,11 +49,44 @@ struct filter_struct { }; typedef struct filter_struct filter; +/* + * A bank of filters that all get applied to the same image. + * + * The point of this is the forward FFT. Applying N filters with N separate + * 'filter' structures transforms the image N times, and every one of those + * transforms is identical. Here the image is transformed once and the result + * is reused, so the cost goes from 2N transforms to N+1. + * + * image_fft has to be kept separate from work_fft because convolve() does its + * complex multiply in place, which would destroy the transform we are trying + * to reuse. + */ +struct filter_bank_struct { + int fft_size; + int image_size; + int n_filters; + int x_size; + int y_size; + + double *fft_vector; + + fftw_plan fft_backward; + fftw_plan fft_forward; + + fftw_complex *image_fft; + fftw_complex *work_fft; + fftw_complex **psf_ffts; +}; +typedef struct filter_bank_struct filterBank; + /* Function Declarations */ void cleanup(filter *); +void cleanupBank(filterBank *); void convolve(filter *, double *, double *); +void convolveBank(filterBank *, double *, double *); void convolveMemo(filter *, double *, double *); filter *initialize(double *, double, int, int, int); +filterBank *initializeBank(double *, int, int, int, int); /* Functions */ @@ -214,6 +247,125 @@ filter *initialize(double *psf, double max_diff, int x_size, int y_size, int est return flt; } + +/* + * cleanupBank() + * + * bank - A pointer to a filterBank structure. + */ +void cleanupBank(filterBank *bank) +{ + int i; + + for(i=0;in_filters;i++){ + fftw_free(bank->psf_ffts[i]); + } + free(bank->psf_ffts); + + fftw_free(bank->fft_vector); + fftw_free(bank->image_fft); + fftw_free(bank->work_fft); + + fftw_destroy_plan(bank->fft_backward); + fftw_destroy_plan(bank->fft_forward); + + free(bank); +} + + +/* + * convolveBank() + * + * Convolve image with every psf in the bank. This is the same arithmetic + * that convolve() does, in the same order, just with the forward transform + * hoisted out of the loop. + * + * bank - A pointer to a filterBank structure. + * image - The image (must be the same size as the psfs the bank was made with). + * results - Pre-allocated storage for n_filters images, in psf order. + */ +void convolveBank(filterBank *bank, double *image, double *results) +{ + int i; + + /* Compute FFT of the image, once. */ + ftmDoubleCopy(image, bank->fft_vector, bank->image_size); + fftw_execute(bank->fft_forward); + + for(i=0;in_filters;i++){ + + /* Multiply by FFT of this psf and compute inverse FFT. */ + ftmComplexMultiply(bank->work_fft, bank->image_fft, bank->psf_ffts[i], bank->fft_size, 0); + fftw_execute(bank->fft_backward); + + /* Copy into the matching slice of result. */ + ftmDoubleCopy(bank->fft_vector, results + i*bank->image_size, bank->image_size); + } +} + + +/* + * initializeBank() + * + * Set things up for FFT convolution with several psfs at once. + * + * psfs - n_filters psfs, contiguous, each (x_size, y_size). + * n_filters - The number of psfs. + * x_size - the size of the psfs in x (slow dimension). + * y_size - the size of the psfs in y (fast dimension). + * estimate - 0/1 to just use an estimated FFT plan. + */ +filterBank *initializeBank(double *psfs, int n_filters, int x_size, int y_size, int estimate) +{ + int i; + double normalization; + filterBank *bank; + + bank = (filterBank *)malloc(sizeof(filterBank)); + + bank->fft_size = x_size * (y_size/2 + 1); + bank->image_size = x_size * y_size; + bank->n_filters = n_filters; + bank->x_size = x_size; + bank->y_size = y_size; + + normalization = 1.0/((double)(x_size * y_size)); + + /* Allocate storage. */ + bank->fft_vector = (double *)fftw_malloc(sizeof(double)*bank->image_size); + bank->image_fft = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*bank->fft_size); + bank->work_fft = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*bank->fft_size); + + bank->psf_ffts = (fftw_complex **)malloc(sizeof(fftw_complex *)*n_filters); + for(i=0;ipsf_ffts[i] = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*bank->fft_size); + } + + /* + * Create FFT plans. This has to happen before anything is written into + * fft_vector, as FFTW_MEASURE overwrites its arrays while planning. + * + * One pair of plans serves the whole bank, rather than a pair per psf. + */ + if (estimate){ + bank->fft_forward = fftw_plan_dft_r2c_2d(x_size, y_size, bank->fft_vector, bank->image_fft, FFTW_ESTIMATE); + bank->fft_backward = fftw_plan_dft_c2r_2d(x_size, y_size, bank->work_fft, bank->fft_vector, FFTW_ESTIMATE); + } + else { + bank->fft_forward = fftw_plan_dft_r2c_2d(x_size, y_size, bank->fft_vector, bank->image_fft, FFTW_MEASURE); + bank->fft_backward = fftw_plan_dft_c2r_2d(x_size, y_size, bank->work_fft, bank->fft_vector, FFTW_MEASURE); + } + + /* Compute FFT of each psf and save. */ + for(i=0;iimage_size, bank->fft_vector, bank->image_size); + fftw_execute(bank->fft_forward); + ftmComplexCopyNormalize(bank->image_fft, bank->psf_ffts[i], normalization, bank->fft_size); + } + + return bank; +} + /* * The MIT License * diff --git a/storm_analysis/sa_library/matched_filter_c.py b/storm_analysis/sa_library/matched_filter_c.py index 06bf5d12..7e3e829e 100644 --- a/storm_analysis/sa_library/matched_filter_c.py +++ b/storm_analysis/sa_library/matched_filter_c.py @@ -33,6 +33,20 @@ m_filter.initialize.restype = ctypes.c_void_p +m_filter.cleanupBank.argtypes = [ctypes.c_void_p] + +m_filter.convolveBank.argtypes = [ctypes.c_void_p, + ndpointer(dtype = numpy.float64), + ndpointer(dtype = numpy.float64)] + +m_filter.initializeBank.argtypes = [ndpointer(dtype = numpy.float64), + ctypes.c_int, + ctypes.c_int, + ctypes.c_int, + ctypes.c_int] + +m_filter.initializeBank.restype = ctypes.c_void_p + class MatchedFilterException(Exception): @@ -88,6 +102,60 @@ def convolve(self, image): return result +class MatchedFilterBank(object): + """ + Convolve one image with several psfs. + + This exists because applying N psfs to the same image with N MatchedFilter + objects computes the same forward FFT of that image N times. Here it is + computed once and reused, so the transform count per call goes from 2N to + N+1. + + The results are the same numbers that N MatchedFilter objects produce, in + psf order. There is no memoization, the caller is expected to be applying + the bank to an image that has changed. + """ + def __init__(self, psfs, fftw_estimate = False): + """ + psfs - A list of psfs, or an (n, x, y) array. All the same shape. + """ + psfs = numpy.asarray(psfs, dtype = numpy.float64) + if (psfs.ndim != 3): + raise MatchedFilterException("psfs must be (n, x, y), got " + str(psfs.shape)) + + self.n_filters = psfs.shape[0] + self.psf_shape = (psfs.shape[1], psfs.shape[2]) + + rc_psfs = numpy.zeros(psfs.shape) + for i in range(self.n_filters): + rc_psfs[i,:,:] = recenterPSF.recenterPSF(psfs[i,:,:]) + rc_psfs = numpy.ascontiguousarray(rc_psfs, dtype = numpy.float64) + + self.mfilter = m_filter.initializeBank(rc_psfs, + self.n_filters, + self.psf_shape[0], + self.psf_shape[1], + int(fftw_estimate)) + + def cleanup(self): + m_filter.cleanupBank(self.mfilter) + self.mfilter = None + + def convolve(self, image): + """ + Returns an (n_filters, x, y) array, one convolution per psf. + """ + if (image.shape[0] != self.psf_shape[0]) or (image.shape[1] != self.psf_shape[1]): + raise MatchedFilterException("Image shape must match psf shape! " + str(image.shape) + " != " + str(self.psf_shape)) + + image = numpy.ascontiguousarray(image, dtype = numpy.float64) + results = numpy.zeros((self.n_filters, self.psf_shape[0], self.psf_shape[1]), + dtype = numpy.float64) + m_filter.convolveBank(self.mfilter, image, results) + + return results + + if (__name__ == "__main__"): import tifffile diff --git a/storm_analysis/test/test_matched_filter.py b/storm_analysis/test/test_matched_filter.py index 260f4ec2..c782912d 100644 --- a/storm_analysis/test/test_matched_filter.py +++ b/storm_analysis/test/test_matched_filter.py @@ -179,6 +179,88 @@ def test_matched_filter5(): flt.cleanup() + +def test_matched_filter_bank1(): + """ + Verify that a bank gives exactly what the same psfs give one at a time. + + This is the property the bank exists to preserve: it only hoists the + forward FFT out of the loop, so every number should be identical, not + merely close. + """ + x_size = 80 + y_size = 90 + n_filters = 4 + + psfs = [] + objects = numpy.zeros((1, 5)) + for i in range(n_filters): + objects[0,:] = [x_size/2, y_size/2, 1.0, 1.0 + 0.4*i, 1.3 + 0.2*i] + psf = dg.drawGaussians((x_size, y_size), objects) + psfs.append(psf/numpy.sum(psf)) + + numpy.random.seed(42) + image = numpy.random.poisson(100, (x_size, y_size)).astype(numpy.float64) + + flts = [matchedFilterC.MatchedFilter(psf, fftw_estimate = True) for psf in psfs] + expected = numpy.array([flt.convolve(image) for flt in flts]) + + bank = matchedFilterC.MatchedFilterBank(psfs, fftw_estimate = True) + result = bank.convolve(image) + + assert (result.shape == (n_filters, x_size, y_size)) + assert numpy.array_equal(expected, result) + + for flt in flts: + flt.cleanup() + bank.cleanup() + + +def test_matched_filter_bank2(): + """ + Verify that a bank of one is just a matched filter, and that mismatched + shapes are rejected rather than read past the end of the array. + """ + x_size = 40 + y_size = 50 + + objects = numpy.zeros((1, 5)) + objects[0,:] = [x_size/2, y_size/2, 1.0, 2.0, 2.0] + psf = dg.drawGaussians((x_size, y_size), objects) + psf = psf/numpy.sum(psf) + + image = numpy.zeros((x_size, y_size)) + image[int(x_size/2), int(y_size/2)] = 3.0 + + flt = matchedFilterC.MatchedFilter(psf, fftw_estimate = True) + bank = matchedFilterC.MatchedFilterBank([psf], fftw_estimate = True) + + assert numpy.array_equal(flt.convolve(image), bank.convolve(image)[0]) + + okay = False + try: + bank.convolve(numpy.zeros((x_size + 2, y_size))) + except matchedFilterC.MatchedFilterException: + okay = True + assert okay, "A wrongly shaped image should raise." + + okay = False + try: + matchedFilterC.MatchedFilterBank(psf) + except matchedFilterC.MatchedFilterException: + okay = True + assert okay, "A 2D psf array should raise." + + flt.cleanup() + bank.cleanup() + + if (__name__ == "__main__"): + test_matched_filter1() + test_matched_filter2() + test_matched_filter3() + test_matched_filter4() test_matched_filter5() + test_matched_filter_bank1() + test_matched_filter_bank2() From 8b98701a43e7cb9d3ddbba2d15b8471aaef3384e Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Sat, 12 Sep 2026 19:09:53 -0400 Subject: [PATCH 2/3] Use a filter bank in the spliner peak finder PeakFinderArbitraryPSF convolves at each of the z values in z_value, and both convolution inputs are the same at every one of them: bg_var for the variance filters, and image - background - fit_peaks for the foreground filters. Each z value was transforming those same two arrays again. With one bank per filter type the two inputs are transformed once per call instead of once per z value, so a run with N z values goes from 4N transforms per peak finding pass to 2N + 2. The foreground subtraction was also being recomputed inside the loop to produce the same array every time, so it is now done once where the convolution happens. Measured on test_spliner.dax with test_spliner_dh.xml, which uses five z values, timing peakFinder() only and alternating the two versions so that neither arm gets a whole block of warm cache: bank 0.176s master 0.235s 1.33x bank 0.170s master 0.237s 1.39x bank 0.172s master 0.291s 1.69x bank 0.219s master 0.291s 1.33x The convolutions themselves go 1.67x at N = 5; peak finding as a whole gains less because findMaxima and the masking are unchanged. That works out to roughly a sixth of steady state analysis time for this configuration. At the default z_value of [0.0] there is one filter per bank and nothing changes, which is the expected result rather than a disappointment - the redundancy only exists when there is more than one z value. Verification. Analyzing with fftw_estimate = 1 gives a byte for byte identical HDF5 before and after, across all 234 datasets. That qualifier matters. With the default fftw_estimate = 0 the output is NOT reproducible, on master or on this branch: two master runs of the same file differ in 2 to 8 datasets at up to 1.1e-15 relative, which is larger than the 5.0e-16 that separates master from this branch. FFTW_MEASURE picks its algorithm from plan time benchmarks, and that choice varies with timing noise. Pinning the planner removes the variation from both arms and leaves them identical, which is the real statement that this change is arithmetically inert. Full suite green at 280 tests. --- storm_analysis/sa_library/fitting.py | 49 +++++++++++++++++----------- 1 file changed, 30 insertions(+), 19 deletions(-) diff --git a/storm_analysis/sa_library/fitting.py b/storm_analysis/sa_library/fitting.py index c30cf27a..72fc7653 100644 --- a/storm_analysis/sa_library/fitting.py +++ b/storm_analysis/sa_library/fitting.py @@ -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 = [] @@ -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): """ @@ -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): """ @@ -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. # @@ -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 From eec2450a371f11ff349833fcf89e03ef36f5b717 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Sat, 12 Sep 2026 19:14:44 -0400 Subject: [PATCH 3/3] Use filter banks in the multiplane peak finder Same redundancy as the spliner peak finder, one dimension over. The filters were stored indexed by z value and then by channel, and the loops ran z on the outside, so each channel's image was transformed once per z value even though the image does not depend on z. The filters are now one bank per channel, each holding that channel's psf at every z value, and each channel's image is convolved once per pass. For N z values and C channels that takes the transform count from 4NC per pass to 2C(N + 1). The loops still accumulate z on the outside and channel on the inside, so the order the per-channel contributions are summed in is unchanged. That matters because those sums are floating point and reordering them would perturb the result. MPPeakFinderDao is left on the same code path but has a single z value, so its banks hold one psf each and it does exactly what it did before. That is expected: like 3D-DAOSTORM, it convolves each image once, and there is no shared transform to hoist. Verified against the multiplane_spline diagnostic, two channels and three z values, with fftw_estimate = 1 so that the FFTW planner cannot vary between runs: all 191 datasets in the output HDF5 are identical before and after. Timing the peak finding step only, alternating the arms: bank 0.539s master 0.699s 1.30x bank 0.635s master 0.684s 1.08x bank 0.627s master 0.676s 1.08x Less than the 1.5x the convolutions themselves get at three z values, because multiplane peak finding also does the affine transforms into the channel 0 frame and the accumulation, neither of which changed. Full suite green at 282 tests. --- storm_analysis/multi_plane/fitting_mp.py | 112 +++++++++++++---------- 1 file changed, 62 insertions(+), 50 deletions(-) diff --git a/storm_analysis/multi_plane/fitting_mp.py b/storm_analysis/multi_plane/fitting_mp.py index 520771a2..abee7070 100644 --- a/storm_analysis/multi_plane/fitting_mp.py +++ b/storm_analysis/multi_plane/fitting_mp.py @@ -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 = [] @@ -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)): @@ -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: @@ -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) @@ -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. @@ -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: @@ -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, @@ -469,14 +486,8 @@ 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: @@ -484,6 +495,13 @@ def setVariances(self, variances): 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) @@ -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: