Abel inversion using Fourier-Hankel Method

I am trying to implement the Abel Inversion of interference microscopy data using the Fourier-Hankel Method. In this method, the Abel Inversion is represented by an inverse Hankel transform following a Fourier transform of the measured data (see [1, 2]), which is the 2D projection of a 3D cylindrically symmetric object.

While performing the fast Fourier transform is easy to do in MATLAB, there seems to be no built-in function to do the Hankel transform. I found a library online [3, 4] which implements the quasi discrete Hankel transform and its inverse. However, I am having troubles using this to recover the original data from a projection.

The following is a minimal example to show my problem. The measured data is a 1D projection (onto the x axis) of a radially symmetric Gaussian function. The reconstruction should be similar to the original, non-projected data. I added a scaling factor of 2 pi / N to the reconstruction, so now the result does not depend on the number of sampling points anymore. However, the amplitude still scales with the outer radius X of the field of view, and also with the width B of the object (see Fig. 1). Additionally, there is a slight difference in shape: Where the original is a Gaussian, the reconstruction seems to be shaped more like a Bessel function, i.e. there is not only a vertical scaling but also a horizontal one (see Fig 2, where the only thing that changed is the value of B).

I am guessing that there is some issue between the definitions of the discrete Fourier and the discrete Hankel transform? Any input would be highly appreciated.

%%Create data for a Gaussian
X = 10;                             % edge of field of view
A = 1;                              % Gaussian amplitude
B=1/sqrt(log(2));                   % Gaussian width
xx = linspace(-X,X,500);            % x axis
f = @(x,z) A*exp(-(x.^2+z.^2)/B^2); % equation for radially symmetric gaussian

dataOriginal = f(xx,0);
dataProjected = integral(@(z) f(xx,z), -X,X, 'ArrayValued', true);

% We could also perform the integration analytically
% fint = @(x) exp(-x.^2/B^2)*sqrt(pi)*B*erf(X/B);
% dataProjected = fint(xx);

figure
hold all
plot(xx, dataOriginal, 'DisplayName', 'original');
plot(xx, dataProjected, 'DisplayName', 'measured projection');

%%Abel inversion
N = floor(length(dataProjected)/2); % Number of sampling positions with positive frequency
dataFFT = real(fft(fftshift(dataProjected))); % The FFT of a symmetric even real function should be real and even.  

% Initialization: retrieve integration kernel I, spatial and frequency
% domain sampling positions r, k and scaling factors R, K
[H,k,r,I,K,R]=dht([],X,N,0); 
% Perform Hankel transform of the fft data
dataReconstructed = idht(dataFFT(1:N)', I, K, R)'./N*2*pi;

plot(r, dataReconstructed, 'DisplayName', 'reconstruction');
legend show
hold off

[1] http://en.wikipedia.org/wiki/Abel_transform#Relationship_to_the_Fourier_and_Hankel_transforms

[2] Montgomery Smith, L., Keefer, D.R., and Sudharsanan, S.I. (1988). Abel inversion using transform techniques. Journal of Quantitative Spectroscopy and Radiative Transfer 39, 367–373.

[3] http://de.mathworks.com/matlabcentral/fileexchange/13371-hankel-transform

[4] http://documents.epfl.ch/users/l/le/leuteneg/www/MATLABToolbox/HankelTransform.html

[x h]

I met the same problem, have you solve it?