Code covered by the BSD License  

Highlights from
Abel Inversion Algorithm

5.0 | 9 ratings Rate this file 76 Downloads (last 30 days) File Size: 4.87 KB File ID: #43639
image thumbnail

Abel Inversion Algorithm



26 Sep 2013 (Updated )

Fourier-based reconstruction of an unknown radial distribution assuming cylindrical symmetry.

| Watch this File

File Information

The reconstruction of the radial density distribution of a cylindrically symmetric object is a common task in different area of physics (e.g. plasma physics). Typically, optical measurements of objects like plasma columns or flames are integrated along the line of sight.
To obtain the underlying distribution from a measured projection, the inverse Abel transform has to be calculated.

This submission provides a Fourier-based algorithm which extracts the radial (2D) distribution from a one-dimensional projection measurement. The algorithm has been proposed and published by G. Pretzler (Z. Naturforsch. 46a, 639 (1991)). Compared to earlier approaches towards Abel inversion, this algorithm is relatively insensitive to noise and errors in the determination of the object's center (see G. Pretzler et al. , Z. Naturforsch. 47a, 955 (1994)).

The fundamental idea is to fit the whole measured profile to a set of cos-expansion-based integrals. (In the conventional approach, in contrast, the radial distribution is obtained by starting at the edges and incrementally iterating towards the center - making it more prone to errors.)

UPDATE: There have been some questions on why the Abel inversion result of generic test data seems to be off by a factor 2. I found a simple explanation: The algorithm assumes that the measurement of line-of-sight-integrated data starts on the center axis of your circular object.
In most experimental situations, however, the measurement takes place along the complete object, effectively doubling the optical path length. Therefore, if you want to analyze experimental data obtained by observing circular (cylindrical) objects, please divide the result F_REC by 2 to yield the right result.

Required Products MATLAB
MATLAB release MATLAB 8.0 (R2012b)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (49)
06 Apr 2015 dj327

dj327 (view profile)

Hi Carsten, the problem of isqcuivefit not running is solved. A classic 1*N matrix or N*1 matrix problem. Although I made sure the input of h is 1*N matrix as requested, dimension and h and myfun inside lsqcurvefit(myfun,x0,hn,h) don't I just simply transposed h

Comment only
02 Apr 2015 dj327

dj327 (view profile)

Thanks for your fast reply, sorry I didn't read through the whole description.

02 Apr 2015 Carsten Killer

@dj327: That depends on whether the line of sight is a full secant line (as in most experimental situations) or just a half secant line, which would generally be sufficient for the problem. The code is written for the latter situation (with hindsight maybe not a good choice, but I do not want to change the submitted code for this minor issue for reasons of consistency).

Please note that this issue is already addressed in the last paragraph in the code description above ("UPDATE: ...").

Comment only
02 Apr 2015 dj327

dj327 (view profile)

I believe you have missed a factor of 2. I need to divide f_rec by 2 to get the expected result.
Could you please check?

Comment only
06 Nov 2014 Carsten Killer

@aneps: Just process every row (or column, depending on the symmetry axis in your system) separatey. For example, you can call the abel_inversion function 100 times for the 100 rows/colmuns, using a simple loop.

Comment only
06 Nov 2014 aneps

aneps (view profile)

My data input (H) is a colormap; a matrix of size 100 x 100. When I tried to use this program to perform inverse Abel transform it shows error : "Error using .*
Matrix dimensions must agree" .

I understand the input data set should be 1 X N matrix. How can I use this program in my case where I have an intensity map of mxn matrix (mine is 100x100)? Thank you.

Comment only
20 Oct 2014 surya Narayan

Whenever i mask the data i have does not start from the axis symmetric center but at a position of some pixels off from the center. How to use abel inversion in this case? Please suggest some references.

Comment only
20 Oct 2014 Carsten Killer

Could you please specify what exactly you mean with "masked" ? (In any case, your data has to be smooth enough to comply with the general constraints of the Abel problem.)

Comment only
18 Oct 2014 surya Narayan

I am using the same code for abel inversion of the axis symmetric data. I have a masked region at the center. Abel inversion is giving me a wrong values if i have a mask. What should i do to make the code work? please send me some references if u have any.

18 Oct 2014 surya Narayan  
12 Oct 2014 Carsten Killer

Hi Laurel,
if the symmetry axis is in the center of your image, the image has to be "cut in half" for Abel processing, i.e. the left and right half are calculated independently. Please remember to use something like fliplr for the left half, so that the first element of the input vector corresponds to the center of the original image.
In your case, an input vector therefore seems to contain only 64 data points. That is not much - I am not sure how well the algorithm will perform in this case.
If are still not sure whether you use the function correctly , you can send a sample image to and I will take a look.
Cheers, Carsten

Comment only
12 Oct 2014 Laurel

Laurel (view profile)

I'm trying to perform an Abel inversion on a set photographs (128x800) of a flame issuing from a jet (to analyze for emissions). The cylindrical symmetry is oriented along the longer (800) axis. Per previous comments, I've been trying run the Abel function one row at a time (running from 0 to 128). However, the produced results do not appear to be accurate. I also tried cutting the image in half and started my iteration from the centerline of the jet. I am unsure if I am using the function in the correct manner.
Thank you

Comment only
14 Aug 2014 Carsten Killer

Recently, the journal (Zeitschrift für Naturforschung) has been digitalized. The paper on Abel inversion can be found on page 639 at

1. Actually, noise is accounted for by using a finite number of cos-expansions. The less expansions you use, the more noise is filtered out (real physical information is of course neglected on the same scale as the noise, so you shouldn't use a too little number of expansion.)

2. You are right, center shift is not included in this code.

Comment only
14 Aug 2014 nudt

nudt (view profile)

The programm is helpful, but several factors such as noise and center shift are not taken into account.
Could you please share reference 1 to me so that I can understand the theoretical background? my email is
Thank you

21 Jul 2014 Carsten Killer

It is not necessary to figure out a function. Just use a row of your image (with the symmetry axis at the first element of the vector) as the first input of the main function ("h").

(the code works numerically (not analytically), e.g. the polynomial sample function in GENERATE_TEST_DATA is also discretized onto a grid and is not treated as an analytical function)

Comment only
18 Jul 2014 Dimas Avila

Thanks for the speedy reply Carsten! So for each row of my image should I figure out a function that characterizes that row and input that?

Comment only
18 Jul 2014 Carsten Killer

@Dimas Avila: The polynomial function is just an example given in the function GENERATE_TEST_DATA. (the result of this function is shown in the sample image at the top of this page). This function is only called if you do not provide input data.

The data vector that you want to perform the Abel inversion on is the first input variable in the main function. Please perform the Abel inversion for each row/column (depending on whether your symmertry axis is vertical or horizontal) of your image separately.

Comment only
17 Jul 2014 Dimas Avila

I was trying to perform an Abel inversion on a set of flame photographs I took to look for emission radicals. I was wondering how I can use this program to perform this. I see there is a function which generates a polynomial distribution function and then takes the abel inversion of that function. I just don't know how to do this for an image.

Comment only
08 Jul 2014 Carsten Killer

@Adnan: I don't see any way in which Abel inversion would be related to your problem. You might want to consider some edge detection algorithm (Matlab offers plenty of them) to find the shape of the person in each frame.

Comment only
08 Jul 2014 Adnan Farooq

@Carsten. In my case i have video of a person who is walking in front of camera (we can see the side view of person) and there is nothing at the background. so all the frames consist of each step he took for walking. so i want to apply Abel inversion projection on each frames (similar like making features vector)

Comment only
07 Jul 2014 Carsten Killer

@Adnan: I am not sure what you mean with "finding changes". Of course, in different measurement images, the Abel inversion results of these images will probably differ as well. But please consider that the Abel inversion is only suited for a special situation, i.e. the investigation of symmetrical objects. Of course you can apply the Abel inversion code individually for each row or column, if your data is some kind of line-of-sight integrated measurement of a symmetrical object/system.
If the axis of symmetry is vertical, you have to consider the rows of your images. Please keep in mind that the axis of symmetry is at the start of your input vector (row) , i.e. the first element of the row vector corresponds to the center of your system r=0 and the last element corresponds the the edge, r=R.

Comment only
07 Jul 2014 Adnan Farooq

@ Carsten and Borek, Thanks for your answers. Actually i want to use Abel transform on sequence of images i have to find the changes in that ... so i want to in that case is it possible to apply Abel inversion directly on 2D row wise and column wise..??

Comment only
07 Jul 2014 Borek

Borek (view profile)

@Adnan Farooq

As Carsten wrote, result is a distribution
f(r), where r is the diameter of a symetric circle. It means, the distribution is same for angle phi from 0 to 2pi. To plot f(r,phi) in 2D, the same parametrization can be used. If you need something more sofisticated, see assymetric abel inversion, which can be used for assymetric signals for example for simple plot of irradiation profile of plasma in tokamaks.

Comment only
07 Jul 2014 Carsten Killer

@Adnan Farooq: Could you maybe pose a more specific question? To summarize, this code takes your input F(y) and calculates the inverse Abel transform, resulting in the corresponding f(r). (using the standard notation, e.g. used in ).

Please note that this code only performs the calculation for one vector, i.e. one row (or column, depending on your systems's geometry) of your image. If you want to transform a complete 2D image, every row (or column...) has to be processed individually.

Comment only
06 Jul 2014 Adnan Farooq

Hello Carsten

Thank you for submitting for code. I am studying about Abel Transform and want to use for 2D image.. i want to know that does this code supports?

10 Jun 2014 Carsten Killer

@surya Narayan:
1) In your case, only half of the data has to be used, i.e. the symmentry axis is at the beginning of your data vector. From the geometrical point of view, both halves (left and right) are independent. For perfect symmetry, they will of course yield identical results. Therefore, the difference between the left and right half isa good indicator for the quality of symmetry in your system.

2) The radius is important for the scaling of the results. For any arbitrary value of R you will get a qualitatively correct result from the Abel inversion (except for very low values of R << 1 , which should be avoided). But, if you need a quantitatively correct result, you have to use the (outer) radius of your object. (e.g. if you investigated some glowing or light absorbing stuff in a glass cylinder which has a radius of 50 mm , just use R=50 and your result will be scaled to mm units)

Comment only
10 Jun 2014 surya Narayan

I have 2-D unwrapped phase values (axis is at the middle of the data)

1. I am entering the axis symmetric data row wise(vertical symmetry is present) as full row. Is this right or should i give only one half of the data?

2. What is the radius that has to be given as input?

can anyone help me out?

23 May 2014 Sina Tebianian  
22 May 2014 Sina Tebianian

Thanks Carsten, Sorry I wrote by mistake "arrow" instead of array. Anyway my data has a physical range that goes from 0.48 to 1. So I understand that you suggest that I subtract the minimum value of the vector from all its components. Once I run the function, I will add the minimum value to the output data again. I figured the normalization once I went through the mfile. Many thanks for your help.

Comment only
22 May 2014 Carsten Killer

Hi Sina,
about the data range: If you use the plot option of my code, both your original data and the calculated radial profile will be plotted normalized to 1 (to make it easer to visually compare the profiles). If you are interested in the correctly scaled data, please use the output (f_rec=abel_inversion(...) ) and plot it manually. (of course, the code uses your original data for the calculations and not the normalized data shown in the figure).

About your data: I am not sure what you mean with an "arrow". However, the boundary value of your data vector (i.e. the last value) should ideally be zero (or something close). This condition is crucial for the algorithm and is imposed by the basic geometry of the situation. If there is some kind of constant background information behind your data, you might want to consider substracting 0.5 in order to create a data vector from 0 to 0.1 (it depends on your specific situation whether this is legitimate).
I hope this helps a bit. Feel free to ask further questions if anything remains unclear.

Comment only
22 May 2014 Sina Tebianian

Hi Carsten,
Thanks for your useful code. I tried to run the program with an arrow of my data but the plots produce data that are completely unexpected. The measured profile that must reflect my data seems to be shifted up. My data range goes from 0.51 to 0.6 but then on the profile I see that the measured profile goes from 0.8 to 1. Also the reconstructed profile is completely different from the measured one. I am wondering if I have not understood part of the code. My vector is basically 1x301 that goes from 0.5 to 0.6 therefore the radial profile should not be very different from these values.
I appreciate your help

Comment only
22 May 2014 Mohammed Ibrahim

@Carsten Killer

I am not able to access reference 1, based on which the program is created. It would be great if you can share it with me.

Comment only
15 May 2014 Carsten Killer

@Mohammed Ibrahim: The "measured distribution" is F(y) and the "reconstructed distribution" is f(r).

(I am sorry for any confusion caused. I thought these denotation for F(y) and f(r) to be practical description for applications in data analysis of physics experiments)

Comment only
15 May 2014 Mohammed Ibrahim

@carsten kiler

when i run the abel inversion for a test case i get measured and reconstructed distribution F(y). How do i get the radial distribution f(r)?

Comment only
15 Feb 2014 Erzberger  
24 Jan 2014 Carsten Killer

By default, lsqcurvefit is not called unless you set the fifth input parameter to 1. If you try the same dataset without using lsqcurvefit (e.g. setting the fifth input to 0), do you also get an error?
Furthermore, do you also get this error when using the lsqcurvefit routine with the built-in test function (i.e. calling "abel_inversion([],[],10,1,1) ")?

Comment only
23 Jan 2014 Akash

Akash (view profile)

Dear Sir
Thanks, Now i am suffering from another kind of problems like

??? Error using ==> lsqcurvefit at 253
Function value and YDATA sizes are

Error in ==> solve_lsq at 17

Error in ==> abel_inversion at 79
A=solve_lsq(h,hn,x0); % solve for
amplitudes A
please help to understand the things.

Thanks in advance

Comment only
22 Jan 2014 Carsten Killer

Hello Akash,
I recently found out that the integral function was introduced in Matlab2012 and is not available in older versions.
However, simply replacing "integral" by "quadgk" should do the job in older Matlab versions.

Comment only
22 Jan 2014 Akash

Akash (view profile)

Hello Sir
I am not able to run this program. Their is an error
" Undefined function or method 'integral' for input arguments of type 'function_handle"
I think integral function is missing. Since i am new to matlab so not 100% sure what is going on. Please help me.
Thanks in advance

Comment only
29 Nov 2013 Carsten Killer

@R: I will try it with Octave as soon as possible. With Matlab, your example works for me: If I do this
and then start the function:
The result is in my case a more or less constant line at y=2 for 0<x<6 (not minding the numerical insecurities now), which is exactly what one would expect.

From the list in Wolfram Mathworld, you might expect the result to be a constant line at y=1 instead of y=2. This difference of the factor 2 is just the consequence of considering the complete half circle or just one quarter segment of the circle (see e.g. in the sketch in )

So, your line at $y approx 2$ should be the right result (or at least the approximation to the analytical result y=2 constant). What do you mean with "a few other messy lines"?

Comment only
28 Nov 2013 R

R (view profile)

Does it work with Octave? I tried a couple of analytic function/inversion pairs and it's not even close. For instance: f(x)=2*sqrt(a^2-x^2); Wolfram Mathworld says the inversion sould be PI_a(x), where PI_a(x)=0<x<a? 1:0. With a=6 and number of cos-expansions=10, on a reconstructed distribution graph, I get a line $x \approx 2$ for 0<=x<=6 and a few other messy lines.

Comment only
28 Nov 2013 Carsten Killer

@aneps: OK, I think now I see what your problem is about. The program will not work with scattered data, i.e. a list of x-y-pairs. You will have to bin your data (map it onto a 2D-grid), e.g. with bilinear interpolation.
After doing this, you will have a 2D-matrix which is some kind of a density plot of your measurement (like taking a picture with a camera).
If y is the axis of symmetry, each row of the matrix can be used as input H to obtain the radial distribution at that particular heigth (y-coordinate).

(note that if the axis of symmetry is located in the center of your matrix and not at x=0, you just have to use "half" of the rows. In an ideal measurement, with a perfect symmetry of your observable, the left and right half of the matrix would of course be identical. But in reality this will probably not be the case.)

The radius R is only important if one wants the result to be scaled to the units suiting your experimental setup. If you just want to check wether the algorith works for you, just set R=10 or any other arbitrary value.

The upper frequency defines how many cos-expansions will be used. (or, more precisely, up to which order the expansion is calculated). A low value for UPF will lead to a relatively rough approximation of your data, but it is on the other hend insensitive to (high-frequency) noise. Choosing higher values for UPF will result in a more and more exact calculation, which accounts for ever smaller nuances in your data.

Comment only
27 Nov 2013 aneps

aneps (view profile)

@Carsten Killer: Thank you for the reply. I obtain the data by projecting the electron wave function on a plane detector (electron velocity map imaging). The axis of symmetry is Y axis (the polarization axis of the laser used). This is the 2D projection of a 3D wave function. The recorded data (points on which the electron collapses for each measurement) is a 2D distribution of points. My program records the coordinates (x,y) of each points. So, I have a big data of x,y values. SO, here how should I give the input arguments? I mean, how should I choose H, Radius R(which is unpredictable in this case), upper frequency etc?

Comment only
27 Nov 2013 Carsten Killer

@aneps: How is your data distributed along the x- and y-direction? If y is the axis of symmetry (e.g. the vertical coordinate in an upright standing cylindrical object), then I assume you want to perform the Abel inversion in x-direction.
Just take each row of your image and compute the Abel inversion for each row (i.e. for each y-value, you can of course do this in a loop) seperately, because the algorithm accepts only 1D-vectors as input.

Concerning the polar coordinates: You can put your measured data (which is in cartesian coordinates, I guess) directly into the algorithm.
However, if the axis of symmetry in your dataset is not along the x- or y-axis, you would have the perform some kind of transformation to produce a dataset that has an "easy" symmetry (i.e. in x- or y-direction).
Feel free to ask for more details, if my answer was not helpful enough :)

Comment only
27 Nov 2013 aneps

aneps (view profile)

I have a data points with X and Y coordinates. Can I use this algorithm to find the abel inversion of the data? Should I convert the data into polar coordinates? Could you please advise me how to use this algorithm to get the abel inversion for my data?

Comment only
23 Oct 2013 James

James (view profile)

I tested this inversion against 3 known analytical abel pairs and found good agreement with this numerical solution. I found that with smooth data with good signal/noise, increasing upf provided a better answer.

23 Oct 2013 James

James (view profile)

I was able to get the code working by replacing lsqcurvefit with a standard non-linear equations solver and simply supplying the residuals. Everything seems to be ok but I still need to through a few theoretical test curves at it. Thanks!

Comment only
22 Oct 2013 Carsten Killer

@James: Thank you for your suggestion, which is quite helpful. I changed that part of the algorithm - now the optimization toolbox is no longer necessary.

Comment only
21 Oct 2013 James

James (view profile)

Would it be possible to use this code without having the optimization toolbox (lsqcurvefit function)? Does it simply find a scaled coefficient for each of the vectors?

Comment only
23 Oct 2013

Instead of using lsqcurvefit (which requires the Optimization Toolbox), the final equation system can now be solved algebraically

18 Feb 2014

Description has been updated. A paragraph regarding the interpretation of the result with respect to the real (physical) units in experimental situations has been added.

Contact us