Incorporate polar onion peeling algorithm into PyAbel

[DanHickstein]

It would be nice to implement one of the inverse Abel transform methods that uses a polar basis set to allow for smooth transforms of noisy data. Polar onion peeling (POP) is one of the most popular methods for this.

There is a nice Matlab implementation of the POP algorithm here. Adi Natan, the author of that code, said that it would be fine if we ported it to python and incorporate it in to PyAbel. He also has the matlab code that generates the basis sets, so we would not be limited to the 1024x1024 basis set that comes with the program.

EDIT: Adi added the code for the basis set generation to the MATLAB site: http://www.mathworks.com/matlabcentral/fileexchange/41064-polar-onion-peeling/content/bld_lut.m

[rth]

That would be great! It would be interesting to see how different implementations behave with noisy data.

[DanHickstein]

Now that #270 has added rBasex, there is probably little motivation to also add pBasex.

[MikhailRyazanov]

Now I need to say here that POP is not really correct. Namely, this procedure:

that generates the projection of an anisotropic ring is valid only if the ring is infinitely narrow (that is, b_r(R) → Abel[δ(R − r)]). In reality, for extended basis functions with different angular dependences, their projections have different radial parts, not just by the (R/r)ⁿ factor.

Here is a comparison of the actual projections and their POP approximations for the triangular basis functions used in rBasex (the POP article is not clear about the shape of their basis functions, but they are unlikely to be narrower than this):

p_{6; 0}(R) here is what they call b₆(R) — projection of the isotropic basis function centered at r = 6, and p_{6; n}(R) are projections of anisotropic basis functions with cosⁿ θ angular dependences.
As can be seen, all POP approximations have noticeable errors, larger for higher angular orders.

However, for larger radii r, the relative widths of the basis functions decrease, so they become effectively narrower, and the POP approximation works better. Perhaps, for real images (with min r ≳ 100 px), the results will be reasonable enough (plus, the authors say that their method has some other problems at low radii).

One potential advantage is that they use Legendre polynomials instead of cosine powers, so the method should be more numerically stable for high angular orders. Although I don't know whether this is really true and whether this higher stability overweighs the increasing approximation errors. If somebody wants to compare, it would be interesting to see the results...

[DanHickstein]

Okay, this is interesting stuff. Is your conclusion is the same as mine: incorporating POP might be interesting, but it's a fairly low priority to incorporate it?

[MikhailRyazanov]

Yes, maybe label this issue as low-priority, but I wouldn't close it.

[adinatan]

Just to add that unlike the paper mentioned, the error that is accumulated closer to the origin at small radii is dealt in the matlab code by restricting the leading beta order that can be used, which created an effective filter for the angle resolution, as well as attenuate the signal the closer it gets to the origin. In addition, this behavior is negligible for radius ~ >100 pixels, so it can be easily circumvented by resampling the domain needed such that the error is restricted to an unimportant neighborhood around the origin.

[MikhailRyazanov]

I guess, if this is ever implemented, it would be nice to have the reference implementation and add any possible improvements as options (maybe enabled by default). In any case, examining results from existing implementations beforehand should help in judging what to do and how. Maybe these deviations are negligible in practice.