Whiten Extremal Value

[garrettwrong]

User reported bug with whitening. Let's confirm and see if we can incorporate any fixes from the MATLAB code.

Hi,

We ran into a strange behavior of the whitening process.
The estimated whitening filter seems to have zero at the center, and so when dividing the filter, you get huge values in the whitened image.
In the MATLAB code I have safeguards against division by small values.
I don’t know if the Python code has implemented them.
Attached is a minimal example demonstrating the problem.
(Rename scripy.txt to script.py, since outlook blocks Python files).
Thanks,
Yoel

Added .txt extension to all files for upload.

example.star.txt
script.py.txt
example.mrcs.txt

[janden]

It's strange that we get a zero value at the center here. In theory, these filer estimates should all be positive. All the same, might be good to add a small epsilon before dividing.

[j-c-c]

Just looking at the noise psd estimation in the AnisotropicNoiseEstimator() below. Line 133 sets the center value of the filter. Looks like if noise_psd_est[mid, mid] and mean_est**2 are equal, the center will be zero (unless I'm misunderstanding this).

ASPIRE-Python/src/aspire/noise/noise.py

Lines 112 to 133 in 1bff8d3

	def estimate_noise_psd(self):
	"""
	:return: The estimated noise variance of the images in the Source used to create this estimator.
	TODO: How's this initial estimate of variance different from the 'estimate' method?
	"""
	# Run estimate using saved parameters
	g2d = grid_2d(self.src.L, indexing="yx", dtype=self.dtype)
	mask = g2d["r"] >= self.bgRadius
	mean_est = 0
	noise_psd_est = np.zeros((self.src.L, self.src.L)).astype(self.src.dtype)
	for i in range(0, self.src.n, self.batchSize):
	images = self.src.images(i, self.batchSize).asnumpy()
	images_masked = images * mask
	_denominator = self.src.n * np.sum(mask)
	mean_est += np.sum(images_masked) / _denominator
	im_masked_f = xp.asnumpy(fft.centered_fft2(xp.asarray(images_masked)))
	noise_psd_est += np.sum(np.abs(im_masked_f**2), axis=0) / _denominator
	mid = self.src.L // 2
	noise_psd_est[mid, mid] -= mean_est**2

At any rate, in our whiten() method we element wise divide by the square root of the noise psd, where as in the matlab code the division only takes place if the psd is larger than some threshold value (around lines 66-70), otherwise the values are set to zero.

If we add some small epsilon before dividing won't that result in a very large value? And would that be problematic?

[janden]

Looks like if noise_psd_est[mid, mid] and mean_est**2 are equal, the center will be zero (unless I'm misunderstanding this).

Right but that is quite unlikely to happen. They are only equal if the mean value of each image is constant (which it may be due to previous preprocessing, i.e. normalize_bg). It's worth checking to see if this is the case.

At any rate, in our whiten() method we element wise divide by the square root of the noise psd, where as in the matlab code the division only takes place if the psd is larger than some threshold value (around lines 66-70), otherwise the values are set to zero.

Yes that corresponds to a pseudoinverse formulation (since the standard inverse does not exist if one of the filter values is zero).

If we add some small epsilon before dividing won't that result in a very large value? And would that be problematic?

Adding an epsilon corresponds to regularizing the inverse (as opposed to replacing it with a pseudoinverse). They're two options with different properties. Note that the pseudoinverse also gives large values if the threshold used is very large (i.e., if your threshold is ε and you have the value ε + 0.001, you'll divide by 1 / (ε + 0.001), which is almost as large as 1 / ε. The regularized option has the advantage that it does not introduce discontinuities in the spectrum while the pseudoinverse may be better suited in other situations.

Before we go further on this. it's important to know why we're getting zeros in the filter estimate. It shouldn't really happen under normal circumstances.

[j-c-c]

@janden Thanks for the clarification. I'll see if I can track down where these zeros are coming from.

[j-c-c]

Ok @janden, I've resumed work on this and I've tracked down where the zeros are coming from. In the estimate_noise_psd() method linked in an above comment we take the centered fft of the masked images. In the case of the script that Yoel provided there is an earlier step of normalize_bg which causes images_masked to have a mean value of zero. In turn, the fft im_masked_f has a center value of essentially zero for each image. This then causes the filter to have a center value close to zero.

[janden]

Ok got it. Thanks for taking a look at this again.

[garrettwrong]

If the only extremal value of concern is actually the center (dc), then we could clamp it with numpy.nan_to_num. That would convert an inf to a large finite number. Assuming the center values are actually 0 the (convolution's) multiply would be 0*large_num. If that works, a reaonably clean implementation would be to use a context for the application of the power and subsequent nan conversion. @j-c-c , let me know if you think that would work; I expect you'll fill me in later. Thanks

[j-c-c]

#1025 added a safeguard against extremal values for power filters with negative power. Closing.

[garrettwrong]

Reopening. It seems something like this may be coming up again after #1125 . I suspect we are running/evaluating something filter related in singles now and that is causing a bit of trouble in practice that we are not covering in unit tests...

Specifically, I found running the gallery/experiments/simulated_abinitio_pipeline.py in single precision fails. The resulting preprocessed images look like they are junk. (Example has been changed to run as doubles for v0.13.0).

Here are the gallery sim images, noise filter and preproc images with doubles:

and the same, as singles... Note the last log line:

(v0.13.0py39) ➜  experiments git:(simgalhotpatch) ✗ python simulated_abinitio_pipeline.py               
2024-10-08 11:30:10,083 INFO [aspire.volume.volume] /Users/gbwright/Library/Caches/ASPIRE-data/emdb_2660.map with dtype float32 loaded as None
2024-10-08 11:30:10,084 INFO [__main__] Original volume map data shape: (1, 360, 360, 360) dtype:float32
2024-10-08 11:30:10,084 INFO [__main__] Initialize CTF filters.
2024-10-08 11:30:10,087 INFO [aspire.source.image] Creating Simulation with 10000 images.
2024-10-08 11:30:10,110 INFO [aspire.source.simulation] Appending CustomNoiseAdder to generation pipeline
2024-10-08 11:30:10,143 INFO [aspire.source.image] Setting max. resolution of source = 32
2024-10-08 11:30:10,241 INFO [aspire.source.image] Caching source images
2024-10-08 11:45:44,057 INFO [__main__] Perform phase flip to input images.
2024-10-08 11:45:44,600 INFO [aspire.source.image] Perform phase flip on source object
2024-10-08 11:45:44,600 INFO [aspire.source.image] Adding Phase Flip Xform to end of generation pipeline
100%|████████████████████████████████████████████████████████████████████████████| 20/20 [00:00<00:00, 28.84it/s]
2024-10-08 11:45:45,478 INFO [aspire.source.image] Whitening source object
2024-10-08 11:45:45,478 INFO [aspire.source.image] Transforming all CTF Filters into Multiplicative Filters
2024-10-08 11:45:45,479 INFO [aspire.source.image] Adding Whitening Filter Xform to end of generation pipeline
2024-10-08 11:45:52,232 WARNING [aspire.operators.filters] PowerFilter setting 1024 extremal filter value(s) to zero.