rBasex, a better pBasex

[MikhailRyazanov]

Here is a new method for analyzing velocity-mapping images, similar to pBasex, but implemented more efficiently.

The main idea is that the Abel transform of f(r) cosⁿ θ is F_n(r) cosⁿ θ, with some other radial part F_n (different for different powers n) but the same angular dependence. And if a suitable basis is chosen for the radial part, the forward transforms f(r) → F_n(r) can be computed analytically for all n. The rest is as usual in all “BASEX”-type methods.

That is, to perform a forward or inverse Abel transform of an image, it is sufficient to extract from it the radial profiles of all cosⁿ θ terms, transform these radial profiles (which is much faster that transforming the whole image) and, if needed, construct the transformed image from the transformed profiles. The transformed radial profiles also directly give various radial distributions needed in VMI studies, thus if the transformed image itself is not needed, the last step can be skipped, making the analysis even faster.

An example of speed comparison with two_point, the fastest of current methods:

PyAbel benchmark run on AMD64 Family 23 Model 17 Stepping 0, AuthenticAMD

time in milliseconds

=== Basis generation ======================================================

Method              n = 513    n = 1025    n = 2049    n = 4097    n = 8193
---------------------------------------------------------------------------
rbasex                 11.7        34.8       116.7       598.9      3248.7
two_point               2.6        13.7        57.7       234.6      1019.7


=== Forward Abel transform ================================================

Method              n = 513    n = 1025    n = 2049    n = 4097    n = 8193
---------------------------------------------------------------------------
rbasex                  2.7        14.6        62.7       302.2      1167.8
rbasex(None)            1.4         6.5        31.7       132.7       516.2


=== Inverse Abel transform ================================================

Method              n = 513    n = 1025    n = 2049    n = 4097    n = 8193
---------------------------------------------------------------------------
rbasex                  2.7        14.7        63.0       302.6      1170.6
rbasex(None)            1.4         6.6        31.8       132.8       516.8
two_point               1.3         9.9        73.7       561.5      4343.1

(here rbasex(None) means that the output image is not constructed)
And, for completeness, the time that must be added to the two_point transform to extract radial distributions from the image:

order = 2, time in milliseconds
=== Full image, cached ======================================================
Method                n = 513    n = 1025    n = 2049    n = 4097    n = 8193
-----------------------------------------------------------------------------
linear
   MIN, none              1.8         7.6        31.0       125.7       498.1

So when only the radial distributions are needed, the rBasex transform is almost free. Otherwise, when the transformed image is also constructed, rbasex is competitive with two_point (and other fast methods) at 1 megapixel (the typical image size) and wins at larger image sizes. The best thing is that the time complexity scales only quadratically with the image radius, in contrast to the cubic scaling in all other methods (except maybe Hansen–Law).

Now let's see how the documentation will be built for this PR...

(More testing, proofreading and fact-checking about pBasex are still neded.)

[MikhailRyazanov]

Addendum: some small changes to the abel.tools.vmi.Distributions class are also included. Some are directly related to rbasex, and some are just usability/consistency improvements.

[stggh]

This is neat! I will examine the method for some VMI images that have an ion-beam (MCP depleted gain) hole.

[DanHickstein]

This is super rad! I will proofread the code more when time allows.

For now, I'm going to close and re-open this PR to see if I can trigger the readthedocs check.

EDIT: Cool, closing and re-opening seems to have the desired effect. Readthedocs build in progress...

[MikhailRyazanov]

@DanHickstein, great! Now the docs are available.

I will proofread the code more when time allows.

Proofread the documentation, please! :–) Grammar, awkward wording and if something isn't clear...

In the meantime I'll think about adding more unit tests (basis, caching, image shape and that higher orders work correctly).

[DanHickstein]

This looks awesome Mikhail! Sorry to take so long to review it. I should have some time over the next few days.

[DanHickstein]

The comparison with the other methods looks nice!

https://gist.github.com/DanHickstein/5874a3bee2dc0eb37cc724049937fcc4

[DanHickstein]

I made it part-way through a review. Will return tomorrow :)

[DanHickstein]

Wow. I am more and more impressed with this work. I was going to ask if you were going to write a paper about it, but it looks like you've already written one! :) https://pyabel--270.org.readthedocs.build/en/270/transform_methods/rbasex-math.html#rbasexmath

[MikhailRyazanov]

Thanks, Dan! I'll wait for you to finish your review (please also make any corrections to the text!) and for Steve's review (@stggh, it's sad that you've got a flu, but hopefully it is less dangerous than this new coronavirus. Get well soon!), and then I'm going to make the edits.

I was going to ask if you were going to write a paper about it, but it looks like you've already written one! :)

Well, this was a practice for a real one, which I hope to publish eventually (about the original method and this faster version). When I get time to work on it...

[stggh]

Amazing and very impressive work @MikhailRyazanov!

The only issue that has appeared in my analysis is a factor of π between the radial intensity yielded by .distr.rIbeta() and that from .angular_integration, no doubt a fault with the latter :-(

im = np.loadtxt('O2-ANU1024.txt.bz2')
imc = abel.tools.center.center_image(im, method='slice')

rbasex = abel.Transform(imc, method='rbasex', angular_integration=True)

rb, rp, rp1 = rbasex.distr.rIbeta()
radial, pes = rbasex.angular_integration

Here is the O₂^- PES, that includes the factor 1/π for the distr intensity (blue line).

The 'extra' intensity is due to the better resolution of the basis decomposition.

The same π scaling factor applies for other species, e.g. C₂H^-:

[stggh]

rbasex sets a new standard for velocity-mapped image analysis. Thanks @MikhailRyazanov for your mammoth effort in implementing, and perfecting this method for PyAbel.

[MikhailRyazanov]

@stggh, we did discuss the π factor in #257. Basically, rbasex is based on Distributions, so it has all the same properties. My I(r) follows the usual convention that 4π ∫ I(r) dr equals the total intensity (the number of events). I don't remember what exactly angular_integration calculates and why. But at least the reconstructed images are consistent with other methods, so you can use either Distributions or angular_integration in any case. :–)

Can you please do C₂H⁻ with reg='pos' and show the low-intensity portion of that plot (say, up to 0.2e5 or 0.1e5)?

[stggh]

#257

Oh, that is right, thanks for the reminder.

you can use either Distributions or angular_integration in any case

Yes, this works very well.

C₂H⁻ with reg='pos'

Note in this examplereg='pos' affects both Distributions and angular_integration. reg produces some nice smoothing.

I will send you the (unpublished) VMI image data via Cloudstor, in case you want to play without the delay of waiting for me to run your nice rbasex analysis.

[MikhailRyazanov]

Thanks for sending the data! Now I see that the negative dip at small radii is due to a dark spot near the image center:

If I mask it, together with the strange contamination in the lower left quadrant, and the strange vertical line (just in case), as described here:

I get more reasonable intensity behavior near the center in both cases — no significant negative intensities beyond the noise level without regularization, and no suspicious hard clamping to zero with the positive regularization:

In any case, it is interesting to see (more obviously in your plots) how the positivity constraints clean up the 300–360 region, such that the peak near 320 becomes apparent. Do you know whether it is real?

[stggh]

Excellent! Thanks for taking a look at the data. The centre line is due to the two CMOS sensors of the camera. The MCP gain depletion results from the ion-beam, in our coaxial spectrometer, which is present during set-up, but switched out, when collecting photoelectrons, although, the gain depletion remains.

You make some good observations, thanks to the (r)basex regularization analysis.

Off-hand I cannot easily say anything about the possibility of a transition near 320 radius. There are many transitions in this region, and it is heavily perturbed. However, I will take a look, in relation to the assignments that have been made. It may take a little time.

[DanHickstein]

Here are a few more comments. I'll continue to work on this when I get a chance tomorrow. still looking awesome!

[DanHickstein]

I made some suggestions for the rBasex documentation. Otherwise, this looks awesome!

[MikhailRyazanov]

Thanks, @DanHickstein! I have made some changes and asked about less trivial ones, please take a look.

[DanHickstein]

Okay, I went through and replied to your questions and comments. Please let me know if I missed any. It looks like there are still some changes that you still need to make.

[MikhailRyazanov]

@DanHickstein, I have more or less implemented your suggestions: restructured “How to use it” to start from Transform; added some words (and a reference) about regularization to “When to use it”; added a reference to onion_peeling with some more details; plus did other, minor corrections.

Most importantly, I have added a check to Transform that raises an error if its origin is used with method='rbasex' in a wrong way. I think, it is better to fail and show a useful message instead of producing unexpected results (with or without a warning).

If everything looks fine, maybe let's merge this? The few remaining issues, like restructuring the documentation, can be addressed later (I'm planning to do some other additions to the method anyway).

[DanHickstein]

@MikhailRyazanov, I had a look at the most recent changes that you've made and they look great!

I think that having Transform raise an error is the right way to go. The error messages are descriptive and should be quite helpful for users.

I think that you can go ahead and merge. Thank you for this fantastic addition!

[MikhailRyazanov]

Thanks!
Merging...

[DanHickstein]

Great work!