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.)
Hi Sung Gyoon Oh,
it is hard to give general advice on a specific problem. The optimum value for upf depends on the kind of noise in your data. I suggest to choose it only as high as you still believe that the small scale structures in the reconstructed profile are real (e.g. if your plasma has no filaments or similar structures one generally expects a relatively smooth emission profile. I also used this code for inversion of plasma emission profiles and typically used upf in the range 20-30).
Regarding the dip at the center: this might be an artifact due to an imperfect estimate for the center position (axis of symmetry) or some slight misalignments in the optical setup (speaking from experience here...)
Hi Carsten, thank you a lot for the codes.
I am trying to reconstruct 2D planar image based on 3D line of sight plasma image which is assumed to have cylinderical symmetry. I tried generating F_rec per each row perpendicular to symmetric axis, and about 150 data points were used for every reconstruction of F_rec. The original h(y) has greatest value at the center and monotonically decay to 0 as it goes to edge, where the profile of h(y) looks like general cos-wave function.However, as I increased the number of upf, which correspond to the number of cos expansions, I found out there seems to be a optimal value around 5~35. Here are the behaviors of F_rec with increasing upf.
At low upf(5~15) the reconstructed profile tend to have profile nearly same as original h(y), where they both look like cos wave. At higher upf(20~30), the profile starts to shows small dip at the center(r=0) with maximum intensity at 0<r_max<R,and at one fixed location. With even higher upf(30~35), the small scale structure start to show up while the overall structure is maintained. For more higher value, noise-looking unsteady peaks are fluctuating on the curves.
I was wondering what is the optimal choices for upf and wanted to know the way I can expect the optimal upf. If you need more information about the explanation above, please send mail to email@example.com. Thank you.
I will look into that.
I am not sure how this problem arises. How do you call the function?
plot_results is the 4th input variable and should be either 1 or 0. If you don't provide it as input when you call the function, the variable is generated in line 44 or 53, depending on whether or not you put in data is first input ("h")
Thank you for posting the code. I can apply directly my data without any extra coding, but when I run the program I get the following error: "Undefined function or variable 'plot_results'." How can I bypass this problem?
Dear Carsten：sorry to bother you again.
my question is to reconstruct 2D chemiluminescence image from a line-of-sight integrated measurement image of a symmetrical flame, hence getting the actual cross profile without accumulating chemiluminescence. I don't know how this code works and where I should change when applying it for my question, could you give me some more specific Instructions to solve my problem. I will be grateful for your reply, my email is firstname.lastname@example.org, thank you very much.
@Carsten Killer:ok,thanks for your reply.
@cooper xi: Yes, the code can be used for this purpose.
If you are unsure how to use the code, you may scroll through the comments below, where you will find some conversations on how to use the code for these applications. If you are still unsure after that, feel free to contact me with specific questions.
Dear Carsten,Thanks for sharing your code.i've got a question concerning the treatment of flame structure, can i use this code to reconstruct 2D images from an measured accumulated flame image? looking forward to your help
i've got a question concerning the normalization of the resulting data.
I figured out, that there is a need of a factor like "length(rec)/3" (rec is the reconstructed data) to fit the initial and the resulting data, but thats not really exact (still about 0.3% mismatch).
Why is there a odd normalization like this necessarry? Or am i mistaken?
Thaks for your time!
Thanks for Crasten,
There is still one problem. In the for loop you wrote where is f_rec and X?
OK, I assume that you have a 1000x1000 image I and the axis of symmetry is vertical in the center of the image, i.e. at x=500. Then, considering the right half of you image (left and right half are independent of each other) you can perform the Abel inversion for line k by:
[f_rec , X] = abel_inversion(I(k,501:end),R,upf);
For the entire image, create a loop through all k and save the result in the k-th row of f
Thanks for your reply. I still have a question. By a for loop you mean the whole program or a specific part. I am not sure where should I define this loop. Would you please tell me where to define this loop. I am really appreciated.
The algorithm only works for one line at a time, so you will have to call it in a loop 1000 times, each time for one distinct row (or column, depending on the direction of your axis of symmetry) of your dataset.
I need to invert my 2D spectrometer images and to do so I can generate 1000* 1000 matrix. But I do not know how to input this matrix. I really appreciate your help.
Thanks a lot! Because my data at the edge is fluctuant, so I cannot so I cannot subtract the background by "h2 = h - h(end)". But I'll try to reduce the effect caused by my data at the edge.
It's really nice to discuss with you!
20 data points is not that much, so the negative result in the center might be an unavoidable result of your limited spatial resolution. Improving the number of data points to 950 by interpolation/fitting does not really add new information, so I am not surprised that this does not change much.
Regarding substracting the background: What I had in my mind does not change the number of data points. If your data is in a 20-element vector "h", just create a new data vector
"h2 = h - h(end)"
which will still contain 20 data points but goes to zero at the edge.
Besides this suggestion, I am afraid that you have to accept the negative results in the center which are probably caused by the spatial resolution of your dataset
Thank you very very much Carsten for the patient and professional respond.
I have 20 data points for one Abel inversion. And If I set Nu as 6,then I got negative reconstructed value in the center but if I Nu is 4, then I won't get negative value.
I also tried to subtract the background, then I only have 16 or 13 data points. but it didn't help too much. I still got negative data even I set Nu as 4.Only that I set Nu smaller than 4.
what's more, I also tried to fit the data first and then subtract the background. then I have about 950 data points. however, even in this case, I can obtain positive value only that I set Nu smaller than 6. Therefore, I think Nu have a significant influence on the value in the center of plasma plume but not the number of my data points. Anyway, I don't need very accurate reconstructed value at the edge because the emission at the edge is naturally low.
in my case, do you still have better suggestions?
How many data points do you have? In any case, Nu should be considerably smaller than the number of data points.
If the data does not fall to zero at the edge you will get wrong reconstruction results at the edge of your plasma since the basic geometry of Abel inversion requires this boundary condition. Whether this error is restricted to the edge region of your data set or whether it might also affect the plasma center depends, again, on the size of your data set and Nu.
As a test, you could manually substract the background light from your dataset, so that the data at the edge becomes zero. Then you can check whether you still get negative emissivities in the center
Morning Carsten, thank you for your reply. I do get a negative value in the center of plasma column once I set the Nu larger than 5. And if I set Nu as 20, then I will get this warning "Matrix is close to singular or badly scaled. Results may be inaccurate." Maybe I just have a few data. Do you think it is a problem that my data doesn't fall to zero at the edge of plasma column?
Dear Wenting He,
the question for the best choice for Nu is indeed not easy to answer. It depends on your problem, e.g. on the number of data points available. While in general a larger Nu allows for more details to be reconstructed you also emphasize noise artifacts and other uncertainties (e.g. a slightly incorrect determination of the radial center of your plasma)
Of course you are right that a negative plasma emissivity is unphysical. Do you still get negative results for larger Nu (e.g. 10,15,20)? Can you rule out that you have other processes influencing your measurement such as re-absorption of light? (such things might also lead to the appearance of negative emissivities).
Hello Carsten, Thanks for sharing your code. It helps a lot. But I still have a question that if I set Nu as 6, I will obtain a negative reconstructed distribution in the center of plasma column. Because my measured profile have a low value in the center of plasma column. I don't think this is correct because the smallest emissivity should be equal or larger than zero. what's the reason for this? PS: if I set Nu as 4 or smaller than 4, the reconstructed distribution won't be negative. But in your explanation, 4 is too low which causes low=pass filtering effect. It looks that the result is highly dependent on the Nu. So how do I know which is the best Nu for my data?
"integral" was introduced in Matlab2012. If you use an earlier version of Matlab, just replace "integral" with "quadgk", which should also work fine.
some functions used in your algorism seems cannot be found, such as 'integral' in compute_expansion.m . Did I should download some other toolbox?
@ot tl: B and L are defined according to equation (8) in the original article by Pretzler, which is in the form A*L=B.
Hence, in my code L is the matrix containing the product of each integral h_n with each other integral h_n. Further, B is defined as the product of the measured data with the integral h_n. On both sides of the equation, we sum up over all positions (index k).
Was this the answer you were hoping for? If not, feel free to ask again.
Hello Carsten, Thanks for sharing your code. I'm trying to go over your code in order to learn how to implement a theoretical mathematical algorithm in MATLAB. However I did not manage to understand the code regarding the solution of the equation system: A*L=B. Can you please explain why did you define the L matrix and B vector the way you did ? ( in particular: what's the meaning of the sums in these definitions ? ) . Thanks a lot !!!
@F. Salehi: Actually the only reason for using this particular basis function is that it was used in the original article by Pretzler (see description), which inspired this code.
You are right, it should principally also work for simpler basis functions. I have not investigated whether the choice of the basis function has any practical impact on the algorithm.
Hi Carsten, Thanks for sharing your scripts here.
I was wondering if there is any specific reason for using (1 - (-1)^n*cos(n*pi*x/R)) as your basis function for expansion instead of simply using cos(n*pi*x/R) ?
This way your basis will be orthogonal.
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 match...so I just simply transposed h
Thanks for your fast reply, sorry I didn't read through the whole description.
@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: ...").
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?
@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.
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.
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.
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.)
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.
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 email@example.com and I will take a look.
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.
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.
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 firstname.lastname@example.org
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)
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?
@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.
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.
@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.
@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)
@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.
@ 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..??
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.
@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 http://en.wikipedia.org/wiki/Abel_transform ).
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.
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?
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)
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?
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.
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.
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
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.
@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)
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)?
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) ")?
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
please help to understand the things.
Thanks in advance
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.
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
@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 https://en.wikipedia.org/wiki/Abel_transform )
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"?
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.
@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.
@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?
@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 :)
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?
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.
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!
@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.
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?
Improvement of the integration scheme (thanks to Benjamin Tadsen). Using the coordinate transformation r=sqrt(t^2+x^2), the full integration limits can be evaluated (in previous versions, an epsilon was substracted to prevent diverging integrals).
I have been made aware of an error in the original article of G. Pretzler which results in missing factor 2. This factor 2 is now included in the code and I have deleted the paragraph in the description above which dealt with this issue.
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.
Instead of using lsqcurvefit (which requires the Optimization Toolbox), the final equation system can now be solved algebraically
Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.