- Email: [email protected]
Abel inversion using transform techniques
J. Quant. Spectrosc.
Radial.
Printed in Great Britain.All
Transfer Vol. 39, No. rights reserved
ABEL
5, Pp. 367-373,1988
0022-4073/88$3.00+ 0.00 Copyright0 1988PergamonPressplc
INVERSION USING TECHNIQUES
TRANSFORM
L. MONTGOMERY SMITH and DENNIS R. KEEFER Center for Laser Applications, University of Tennessee Space Institute, Tullahoma, TN 37388, U.S.A. and S. I. SUDHARSANAN Department of Electrical Engineering, University of Arizona, Tucson, AZ 85721, U.S.A. (Received 8 September
1987)
method is presented for calculating the reconstruction of a circularly symmetric, two-dimensional function from its projection, a relation known as the Abel inversion. This technique differs from earlier methods by using integral transforms for its implementation. The frequency-domain analysis allows for experimentally obtained data, which are often noisy and off-center, to be procesed in a systematic manner. The formulation of the Abel inversion in terms of transforms, filtering of the noise, and estimate of the off-center shift are discussed. Sample calculations of simulated noisy data and the application of the method to an image of a laser-sustained plasma are presented. Abstract-A
1. INTRODUCTION Optical emission or transmission is often used to obtain experimental diagnostic measurements of axially symmetric sources such as plasmas, flames, and rocket and jet engine plumes. The line-of-sight projection of the emission or transmission coefficients onto a one-dimensional axis can be used to determine their two-dimensional radial structure. As is the case in most inverse problems, however, the conditions are ill-posed, and so considerable difficulty can arise when processing real experimentally acquired data that has inherent uncertainties. The reconstruction of a circularly symmetric two-dimensional function from its projection onto an axis is known as Abel inversion or inverse Abel transformation of the projection. The measured intensity, I(X), is given in terms of emission coefficients, c(r), through the Abel transform’ Z(x) = 2
cD K(T) dr (1)
s X J_’ where x is measured symmetric transform,
the displacement of the intensity profile and r is the radial distance in the source. The intensity I(x) is the one-dimensional projection of the two-dimensional, circularly function having c(r) as a radial slice. The inversion integral, or the inverse Abel is given by c(r)=
-;
m (dlldx) 1 I p7
dx .
(2)
In practice, application of Eq. (2) is made difficult because of the singularity in the integral at the lower limit and because the derivative of the projection tends to enhance greatly the noise-corrupting of the data. Furthermore, the intensity is usually not available as a continuous function and so is known only at discrete sample points. Several approaches, based upon geometrical techniques or numerical methods and using polynomial fits, have been employed to perform the inversion. Nestor and Olsen’ transformed the variables according to r2 = u and x2 = u so that the inversion integral can be approximated by a simpler sum. Bockasten3 fitted third-degree polynomials to the data points and approximated the integral by a sum. However, these methods require prior smoothing of the data and are not considered complete in themselves. 367
L.
368
MONTGOMERY SMITH et al
Later, least-squares curve-fitting methods were employed and were found to yield better results than the exact fit methods when applied to noisy data. Freeman and Katz4 used a single polynomial curve to fit the data. Fourth-order polynomials were found to give the best results among trials using up to twelfth-order polynomials. Cremers and Birkebak’ compared several inversion techniques and showed that the least-squares curve fitting techniques are more favorable than the exact fit methods. Short descriptions of several other techniques and comparisons of these various techniques can be found in Ref. 5. Shelby6 divided the data into several intervals and used a least-squares polynomial fit technique in each interval to smooth the scattered data. The inversion was then performed analytically and summed over the intervals to obtain emission coefficients. Malonado et al’ expanded t(r) in a series of orthogonal polynomials and derived a method to find the expansion coefficient from the intensity data. All of these methods have drawbacks. The singularity in the lower limit of the integral causes problems for the numerical methods, but these are avoided by using the analytical methods. The smoothing techniques used are essentially a kind of lowpass filtering having undetermined filter characteristics. When the inversion is performed, the spectral characteristics of the noise, of the desired signal and of the smoothing algorithm are not considered. Therefore, the possible problems of loss of information and distortion are neglected. Use of the Abel integral implies an assumption that the input data will be symmetric, but the measured intensity data often exhibit some degree of asymmetry, and the exact axis of symmetry is not usually known a priori. Determination of the axis of symmetry is often chosen using ad hoc methods. The smoothing techniques consume a large amount of computer time and the error-propagation calculations are tedious6 We present a method based on integral transforms that removes many of the difficulties and uncertainties associated with these earlier techniques and, by using the fast Fourier transform (FFT) algorithm to implement the procedure, greatly reduces the computation time. This technique is substantially different from earlier methods in the use of transform techniques and frequencydomain analysis. Several principles of digital signal processing and spectral analysis are employed to solve the problems associated with the processing of actual, noise-corrupted, asymmetric data. In Sec. 2 the Abel inversion integral is reformulated in terms of the Fourier and Hankel transforms. This reformulation is then used to handle many of the ill-posed conditions encountered when analyzing actual noisy data, and the methods for doing this are outlined in Sec. 3. Numerical experiments demonstrating the validity and applicability of the inversion scheme are presented in Sec. 4, together with a recent application of the method to the image of a laser-sustained plasma. 2.
REFORMULATION
OF
THE
ABEL
INVERSION
The mathematical basis for the Abel inversion technique presented here is the reformulation of the Abel transform equations (1) and (2) in terms of the Fourier and Hankel integral transforms. This formulation is discussed in Bracewell and can also be derived as a special case of the projection-slice theorem, a fundamental relation in the field of computed tomography. For completeness, a brief derivation is presented here. If the substitution Y= 7 x + y is made in Eq. (1) the projection can be written Z(x) = The one-dimensional
Fourier transform 9{Z(x))
If the variables of integration shown that
m c(Jm) s px of Eq. (3) is
dy.
(3)
- i27cxq) dx dy. 3o a: E(,/m)exp( s -cc s -3c are now changed from Cartesian to polar coordinates,
=
P{Z(x)} = 27r ccrc(r)J,,(2nrq) s0
dr,
(4)
it can be
(5)
where J,( .) is the zero-order Bessel function of the first kind. The right-hand side of Eq. (5) is the zero-order Hankel transform of e(r), whose inverse transform is identical in form to the forward
Abel inversion using transform techniques
369
transform. The emission distribution of the plasma can thus be recovered by taking the inverse Hankel transform of the Fourier transform of the projected intensity as follows: c(r) = 2x JornqJ,,(2nrq) [Tm Z(x)exp(-i2zxq)
dx dq.
(6)
From a computational point of view, the inversion formula given in Eq. (6) has several advantages over the Abel inversion integral of Eq. (2). First, we avoid the difficulty associated with the singularity at the lower limit of integration. Second, following the Fourier transform of Z(x), filters may be applied directly in the frequency domain to reduce the noise and thus smooth the data in a known, systematic manner. Finally, Eq. (6) can be numerically approximated with discrete fast transform algorithms available in software or vector-array processing hardware to decrease computation time over techniques used previously. Equation (6) shows an important property of the transform-based formulation of the Abel inversion. If Z(x) is an error-free projection of a real radially symmetric function that has an axis of symmetry that projects onto the origin of the x-coordinate axis, as was implicitly assumed in the derivation, then Z(x) is real and even. Its Fourier transform is thus also real and even, and so the inverse Hankel transform yields a real e(r). If, however, the projection is such that the axis of symmetry is shifted from the x-axis origin and corrupted with some uncertainty called noise, as is the case in any experimental data-acquisition system, then the resulting inversion will have an imaginary component, which is a physical impossibility. It is thus necessary to process any experimentally-acquired data in such a manner so as to ensure that it is symmetric about the x-axis origin when employing Eq. (6) to perform the Abel inversion. This procedure involves determining the axis of symmetry in the presence of noise and eliminating the odd component of the noise from the data. This technique will be discussed in Sec. 3. The computational technique employed for Eq. (6) was to use a fast Fourier transform algorithm to approximate the inner integral by a rectangular rule integration. The inverse Hankel transform was then performed by the method of Candel,’ modified as given in Ref. 10 to increase computational efficiency further by eliminating redundant additions in the algorithm. Numerical experience has shown that this method is quite fast, and results in close agreement between known closed-form Abel inversions and numerically calculated inverted functions. 3.
APPLICATION
TO
NOISY
DATA
The problems associated with the application of the integral transform formulation of Eq. (6) to actual experimentally-acquired data are: (a) how to reduce the effects of noise corrupting the data, (b) how to determine the axis of symmetry, and (c) how to symmetrize the data further about the determined axis. The transform technique we developed allows for all three of these problems to be addressed in the frequency domain following the one-dimensional Fourier transform of the projected intensity data. In this section, the mathematical theory behind this method is presented. The projected intensity measurement incident on the detector is assumed to be shifted from the origin of the detector, corrupted with additive random noise, and sampled at a fixed spatial interval. The resulting signal is thus modeled as s(m)=Z(m
+Am)+n(m),
(7)
where m is the discrete sample number, Am is the shift from the origin of the detector, and n(m) denotes the sampled values of a random noise process. The discrete Fourier transform of Eq. (7) is given by S(k) = G(k)exp(i2nkAm/M)
+ N(k),
(8)
where S(k) and N(k) denote the Fourier transforms of s(m) and n(m), respectively, as a function of the discrete frequency variable k. G(k) is the Fourier transform of Z(m), the unshifted and therefore symmetric intensity projection. The total number of points in the discrete Fourier transform is M. Removal of a significant portion of the noise power is possible by multiplication with the transfer function of a chosen filter to the frequency-domain data of Eq. (8). Because of the baseband nature
310
L. MONTCXXERY SMITH et al
of plasma data, a lowpass filter is suitable for removing the high frequency components of the noise while preserving the signal information contained in the lower frequencies. Filter design can be accomplished through standard optimization programs for desired passband-stopband specifications. Following the initial filtering of the noise, the next step is to determine the shift Am of the data from the origin of the detector. To accomplish this, attention is restricted to the lower frequency components contained within the passband of the previously applied filter where Eq. (8) is still valid. Assuming the noise is stationary, Eq. (8) can be written G(k) = S(k)exp( - i2nkAm/M)
- N(k).
(9)
The function G(k) is the Fourier transform of the real, even function Z(m) and thus is also real and even. Equation (9) can then be broken up into its real and imaginary components as G(k) = S,(k)cos(2~~kAm/M) 0 = S,(k)cos(27ckAm/M)
+ S,(k)sin(2~kAm/M) - S,(k)sin(2nkAm/M)
- NR(k),
(lOa)
- N,(k),
(lob)
where the subscripts R and I refer to the real and imaginary components of each function. Equation (lob) shows that the imaginary portion of the noise spectrum can be written explicitly as a function of the measurement data and the shift of the data from the origin. The noise was further assumed to be an uncorrelated, bandlimited Gaussian random process. Its spectral components can thus be modeled as independent Gaussian random variables,” each with probability density function
~dN,(k))=~ -~N:(k) 202 1 fi exp [
(11)
Substituting Eq. (lob) into Eq. (11) and taking the ensemble probability density functidn over all values of k in the passband of the filter as the product of the individual probability functions yields the conditional probability density of S(k) given Am as
AWlAm)
=
’
(JSco)K
- &
Izd [S,(k)cos(2nkAm/M)
- s,(k)sin(2~kAmiM)12},
exp
(12)
where K is the number of spectral components contained within the passband. Maximizing this probability density function with respect to Am provides a maximum likelihood estimator for the shift. The resulting transcendental equation for the estimate, Am, is K-l C
k{2SR(k)S,(k)cos(4nkAh/M)
- [S;(k) - S:(k)]sin(4xkAfi/M)}*
= 0.
(13)
k=O
This equation can be solved by iterative techniques to determine the estimate for the shift. Once the shift has been determined, the spectrum of the data is multiplied by the phase term exp( - 27rrkAriz/M) to shift the estimated location of the axis of symmetry of the data back to the origin. Following this multiplication, the remaining imaginary components are by assumption necessarily due to noise, and so can be set equal to zero. This process completes the symmetrizing of the data, with the results that the inverse Hankel transform yields a real, even, and smoothed Abel inversion of the data. As a final step in the Abel inversion by means of transform methods, an error analysis based on the statistical properties of the input noise and the effective bandwidth of the applied filter can be made. It can be shown that for wide-sense-stationary input white noise with variance o* filtered with a bandwidth B and inverted via the Abel inversion integral in Eq. (2), the variance of the output noise is approximately’2 var[no(r)] N a2B/2r,
(14)
for sufficiently large r. If the imaginary portion of the noise is eliminated, as described previously, the output noise power is effectively halved and the factor in the denominator in Eq. (14) is 4 for the case of Abel inversion performed by the method presented here.
Abel inversion using transform techniques 4.
NUMERICAL
AND
371
EXPERIMENTAL
RESULTS
A number of numerical experiments have been carried out using both noise-free and noisy data to evaluate the performance of the integral transform method implemented using discrete transforms.‘O These numerical experiments using noise-free Gaussian and semi-elliptic data established the accuracy of the modified Candel method for the Hankel transform. Numerical experiments using noisy data were performed to establish the validity of the method used to symmetrize the data and to evaluate the effectiveness of the optimal filters used to smooth the data. A noisy trial function having a known inversion was generated by sampling a Gaussian function defined as Z,(X) = exp( -x*/400). The Abel inversion of Z,(x) is theoretically
(15)
given by’
6, (r) = (1/20&)exp(
- r*/400).
(16)
A zero mean sequence of random numbers having a Gaussian density function with a variance of 0.001 was generated and added to the sampled sequence of Eq. (15), which was then shifted by two integer sample points. Also, the sequence was zero padded so that the filtering by discrete Fourier transforms produces a linear convolution filtering. This trial sequence is shown in Fig. 1.
-.,ij -100
-
50
0
50
1
0
x
Fig. 1. Input noisy and off-center trial sequence I,(X).
-0.02
5 0
20
40
60
60
100
Fig. 2. Result of the Abel inversion for the sequence shown in Fig. I without filtering or symmetrizing. Circles mark the discrete calculated values connected by dotted lines while the theoretical inversion is shown by the solid line.
312
L. MONTGOMERY SMITHet al 0.04 r
I 20
-0.02 I 0
I 40
I 60
I 60
1 100
r Fig. 3. Result of the Abel inversion of the sequence shown in Fig. I with filtering and symmetrizing of the data. Circles mark the discrete calculated values connected by dotted lines, while the theoretical inversion is shown by the solid line.
The result obtained when this trial sequence was inverted without filtering or symmetrizing is shown in Fig. 2. The same trial function was then inverted using the symmetrizing procedure and an optimal filter having a passband from 0.0 to 0.0525 and a stopband from 0.0725 to 0.5. The result from this inversion is shown in Fig. 3. A principal advantage of this integral transform method is its computational efficiency. It has been applied to the inversions of a large number of images obtained from experiments with laser-sustained plasmas. I3 Each of these digitally sampled images consisted of 512 radial scans of the plasma and each radial scan consisted of 240 sample points. The inversion of 512 sets of 240-point data consumed about 8 min of CPU time when processed by the integral transform method, in contrast to a 20 h CPU time required for a curve-fit inversion using the method described in Ref. 6. The data analysis was carried out using a Masscomp MC500 Micro Supercomputer. The symmetrizing procedure generally consumed about 5-7 iterations for each set of data. Contour plots of the intensity and the emission coefficient for one of these images are given in Figs. 4 and 5, respectively. This particular image illustrates the performance of the method for intensity profiles which exhibit some asymmetry and which result in inversion profiles that are peaked on-axis (e.g. 49 mm), flat-topped (e.g. 45 mm) and having off-axis peaks (e.g. 43 mm).
3.0 2.4 1.8
E 2
1.2
0) 0.6 ,o i
2.4
t
Axial
Fig. 4. Contour
distance
(mm)
plot of a digitally sampled experimental intensity distribution laser-sustained plasma.
obtained
from a
373
Abel inversion using transform techniques
2.4 1.6 z 5
1.2
?j 0.6 ; .b! 0.0 v _ 0.6 0 G 1.2 B 1.6
3.0 40
43
46 Axial
distance
49
52
(mm)
Fig. 5. Contour plot of the plasma emission coefficient obtained by the Abel invesion of Fig. 4 using transform techniques.
5. CONCLUDING
REMARKS
A new integral transform method has been developed for the Abel inversion of measured intensity data in the presence of noise. The integral transform method has a number of advantages compared with existing curve-fit techniques. The curve-fitting methods do not allow for consideration of the spectral characteristics of the noise and the desired information or address the problem of symmetrizing the experimental data. The Abel inversion or the inverse Abel transform is obtained by performing a Fourier transform followed by the inverse Hankel transform. The efficiency of the well known FFT is exploited in the evaluation of both the Fourier and inverse Hankel transforms. This fast technique is made applicable to noisy and off-center experimental data by using frequency-domain filtering and a maximum likelihood estimator derived from the assumed noise characteristics. REFERENCES 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13.
H. R. Griem, Plasma Spectroscopy, McGraw-Hill, New York, NY (1964). 0. H. Nestor and H. N. Olsen, SIAM Rev. 2 2, 200 (1960). K. Bockasten, JOSA 51, 943 (1961). M. J. Freeman and D. Katz, JOSA 53, 1172 (1963). C. J. Cremers and R. C. Birkebak, Appl. Opt. 5, 1057 (1966). R. T. Shelby, Masters Thesis, Department of Mathematics, University of Tennessee, Knoxville, TN (1976). C. D. Maldonado, A. P. Caron, and H. N. Olsen, JOSA 55, 1247 (1965). R. N. Bracewell, The Fourier Transform and Its Applications, McGraw-Hill, New York, NY (1965). S. M. Candel, Comp. Phys. Commun. 23, 343 (1981). S. I. Sudharsanan, Masters Thesis, Department of Electrical Engineering, University of Tennessee, Knoxville, TN (1986). K. S. Shanmugam, Digital und Analog Communications Systems, Wiley, New York, NY (1979). L. M. Smith, IEEE Trans. Inform. Theory, to be published. R. D. Welle, D. R. Keefer, and C. Peters, AIAA Paper 86-1077, “Energy Conversion Efficiency in High-Flow Laser-Sustained Argon Plasmas,” Atlanta (1986).
Radial.
Printed in Great Britain.All
Transfer Vol. 39, No. rights reserved
ABEL
5, Pp. 367-373,1988
0022-4073/88$3.00+ 0.00 Copyright0 1988PergamonPressplc
INVERSION USING TECHNIQUES
TRANSFORM
L. MONTGOMERY SMITH and DENNIS R. KEEFER Center for Laser Applications, University of Tennessee Space Institute, Tullahoma, TN 37388, U.S.A. and S. I. SUDHARSANAN Department of Electrical Engineering, University of Arizona, Tucson, AZ 85721, U.S.A. (Received 8 September
1987)
method is presented for calculating the reconstruction of a circularly symmetric, two-dimensional function from its projection, a relation known as the Abel inversion. This technique differs from earlier methods by using integral transforms for its implementation. The frequency-domain analysis allows for experimentally obtained data, which are often noisy and off-center, to be procesed in a systematic manner. The formulation of the Abel inversion in terms of transforms, filtering of the noise, and estimate of the off-center shift are discussed. Sample calculations of simulated noisy data and the application of the method to an image of a laser-sustained plasma are presented. Abstract-A
1. INTRODUCTION Optical emission or transmission is often used to obtain experimental diagnostic measurements of axially symmetric sources such as plasmas, flames, and rocket and jet engine plumes. The line-of-sight projection of the emission or transmission coefficients onto a one-dimensional axis can be used to determine their two-dimensional radial structure. As is the case in most inverse problems, however, the conditions are ill-posed, and so considerable difficulty can arise when processing real experimentally acquired data that has inherent uncertainties. The reconstruction of a circularly symmetric two-dimensional function from its projection onto an axis is known as Abel inversion or inverse Abel transformation of the projection. The measured intensity, I(X), is given in terms of emission coefficients, c(r), through the Abel transform’ Z(x) = 2
cD K(T) dr (1)
s X J_’ where x is measured symmetric transform,
the displacement of the intensity profile and r is the radial distance in the source. The intensity I(x) is the one-dimensional projection of the two-dimensional, circularly function having c(r) as a radial slice. The inversion integral, or the inverse Abel is given by c(r)=
-;
m (dlldx) 1 I p7
dx .
(2)
In practice, application of Eq. (2) is made difficult because of the singularity in the integral at the lower limit and because the derivative of the projection tends to enhance greatly the noise-corrupting of the data. Furthermore, the intensity is usually not available as a continuous function and so is known only at discrete sample points. Several approaches, based upon geometrical techniques or numerical methods and using polynomial fits, have been employed to perform the inversion. Nestor and Olsen’ transformed the variables according to r2 = u and x2 = u so that the inversion integral can be approximated by a simpler sum. Bockasten3 fitted third-degree polynomials to the data points and approximated the integral by a sum. However, these methods require prior smoothing of the data and are not considered complete in themselves. 367
L.
368
MONTGOMERY SMITH et al
Later, least-squares curve-fitting methods were employed and were found to yield better results than the exact fit methods when applied to noisy data. Freeman and Katz4 used a single polynomial curve to fit the data. Fourth-order polynomials were found to give the best results among trials using up to twelfth-order polynomials. Cremers and Birkebak’ compared several inversion techniques and showed that the least-squares curve fitting techniques are more favorable than the exact fit methods. Short descriptions of several other techniques and comparisons of these various techniques can be found in Ref. 5. Shelby6 divided the data into several intervals and used a least-squares polynomial fit technique in each interval to smooth the scattered data. The inversion was then performed analytically and summed over the intervals to obtain emission coefficients. Malonado et al’ expanded t(r) in a series of orthogonal polynomials and derived a method to find the expansion coefficient from the intensity data. All of these methods have drawbacks. The singularity in the lower limit of the integral causes problems for the numerical methods, but these are avoided by using the analytical methods. The smoothing techniques used are essentially a kind of lowpass filtering having undetermined filter characteristics. When the inversion is performed, the spectral characteristics of the noise, of the desired signal and of the smoothing algorithm are not considered. Therefore, the possible problems of loss of information and distortion are neglected. Use of the Abel integral implies an assumption that the input data will be symmetric, but the measured intensity data often exhibit some degree of asymmetry, and the exact axis of symmetry is not usually known a priori. Determination of the axis of symmetry is often chosen using ad hoc methods. The smoothing techniques consume a large amount of computer time and the error-propagation calculations are tedious6 We present a method based on integral transforms that removes many of the difficulties and uncertainties associated with these earlier techniques and, by using the fast Fourier transform (FFT) algorithm to implement the procedure, greatly reduces the computation time. This technique is substantially different from earlier methods in the use of transform techniques and frequencydomain analysis. Several principles of digital signal processing and spectral analysis are employed to solve the problems associated with the processing of actual, noise-corrupted, asymmetric data. In Sec. 2 the Abel inversion integral is reformulated in terms of the Fourier and Hankel transforms. This reformulation is then used to handle many of the ill-posed conditions encountered when analyzing actual noisy data, and the methods for doing this are outlined in Sec. 3. Numerical experiments demonstrating the validity and applicability of the inversion scheme are presented in Sec. 4, together with a recent application of the method to the image of a laser-sustained plasma. 2.
REFORMULATION
OF
THE
ABEL
INVERSION
The mathematical basis for the Abel inversion technique presented here is the reformulation of the Abel transform equations (1) and (2) in terms of the Fourier and Hankel integral transforms. This formulation is discussed in Bracewell and can also be derived as a special case of the projection-slice theorem, a fundamental relation in the field of computed tomography. For completeness, a brief derivation is presented here. If the substitution Y= 7 x + y is made in Eq. (1) the projection can be written Z(x) = The one-dimensional
Fourier transform 9{Z(x))
If the variables of integration shown that
m c(Jm) s px of Eq. (3) is
dy.
(3)
- i27cxq) dx dy. 3o a: E(,/m)exp( s -cc s -3c are now changed from Cartesian to polar coordinates,
=
P{Z(x)} = 27r ccrc(r)J,,(2nrq) s0
dr,
(4)
it can be
(5)
where J,( .) is the zero-order Bessel function of the first kind. The right-hand side of Eq. (5) is the zero-order Hankel transform of e(r), whose inverse transform is identical in form to the forward
Abel inversion using transform techniques
369
transform. The emission distribution of the plasma can thus be recovered by taking the inverse Hankel transform of the Fourier transform of the projected intensity as follows: c(r) = 2x JornqJ,,(2nrq) [Tm Z(x)exp(-i2zxq)
dx dq.
(6)
From a computational point of view, the inversion formula given in Eq. (6) has several advantages over the Abel inversion integral of Eq. (2). First, we avoid the difficulty associated with the singularity at the lower limit of integration. Second, following the Fourier transform of Z(x), filters may be applied directly in the frequency domain to reduce the noise and thus smooth the data in a known, systematic manner. Finally, Eq. (6) can be numerically approximated with discrete fast transform algorithms available in software or vector-array processing hardware to decrease computation time over techniques used previously. Equation (6) shows an important property of the transform-based formulation of the Abel inversion. If Z(x) is an error-free projection of a real radially symmetric function that has an axis of symmetry that projects onto the origin of the x-coordinate axis, as was implicitly assumed in the derivation, then Z(x) is real and even. Its Fourier transform is thus also real and even, and so the inverse Hankel transform yields a real e(r). If, however, the projection is such that the axis of symmetry is shifted from the x-axis origin and corrupted with some uncertainty called noise, as is the case in any experimental data-acquisition system, then the resulting inversion will have an imaginary component, which is a physical impossibility. It is thus necessary to process any experimentally-acquired data in such a manner so as to ensure that it is symmetric about the x-axis origin when employing Eq. (6) to perform the Abel inversion. This procedure involves determining the axis of symmetry in the presence of noise and eliminating the odd component of the noise from the data. This technique will be discussed in Sec. 3. The computational technique employed for Eq. (6) was to use a fast Fourier transform algorithm to approximate the inner integral by a rectangular rule integration. The inverse Hankel transform was then performed by the method of Candel,’ modified as given in Ref. 10 to increase computational efficiency further by eliminating redundant additions in the algorithm. Numerical experience has shown that this method is quite fast, and results in close agreement between known closed-form Abel inversions and numerically calculated inverted functions. 3.
APPLICATION
TO
NOISY
DATA
The problems associated with the application of the integral transform formulation of Eq. (6) to actual experimentally-acquired data are: (a) how to reduce the effects of noise corrupting the data, (b) how to determine the axis of symmetry, and (c) how to symmetrize the data further about the determined axis. The transform technique we developed allows for all three of these problems to be addressed in the frequency domain following the one-dimensional Fourier transform of the projected intensity data. In this section, the mathematical theory behind this method is presented. The projected intensity measurement incident on the detector is assumed to be shifted from the origin of the detector, corrupted with additive random noise, and sampled at a fixed spatial interval. The resulting signal is thus modeled as s(m)=Z(m
+Am)+n(m),
(7)
where m is the discrete sample number, Am is the shift from the origin of the detector, and n(m) denotes the sampled values of a random noise process. The discrete Fourier transform of Eq. (7) is given by S(k) = G(k)exp(i2nkAm/M)
+ N(k),
(8)
where S(k) and N(k) denote the Fourier transforms of s(m) and n(m), respectively, as a function of the discrete frequency variable k. G(k) is the Fourier transform of Z(m), the unshifted and therefore symmetric intensity projection. The total number of points in the discrete Fourier transform is M. Removal of a significant portion of the noise power is possible by multiplication with the transfer function of a chosen filter to the frequency-domain data of Eq. (8). Because of the baseband nature
310
L. MONTCXXERY SMITH et al
of plasma data, a lowpass filter is suitable for removing the high frequency components of the noise while preserving the signal information contained in the lower frequencies. Filter design can be accomplished through standard optimization programs for desired passband-stopband specifications. Following the initial filtering of the noise, the next step is to determine the shift Am of the data from the origin of the detector. To accomplish this, attention is restricted to the lower frequency components contained within the passband of the previously applied filter where Eq. (8) is still valid. Assuming the noise is stationary, Eq. (8) can be written G(k) = S(k)exp( - i2nkAm/M)
- N(k).
(9)
The function G(k) is the Fourier transform of the real, even function Z(m) and thus is also real and even. Equation (9) can then be broken up into its real and imaginary components as G(k) = S,(k)cos(2~~kAm/M) 0 = S,(k)cos(27ckAm/M)
+ S,(k)sin(2~kAm/M) - S,(k)sin(2nkAm/M)
- NR(k),
(lOa)
- N,(k),
(lob)
where the subscripts R and I refer to the real and imaginary components of each function. Equation (lob) shows that the imaginary portion of the noise spectrum can be written explicitly as a function of the measurement data and the shift of the data from the origin. The noise was further assumed to be an uncorrelated, bandlimited Gaussian random process. Its spectral components can thus be modeled as independent Gaussian random variables,” each with probability density function
~dN,(k))=~ -~N:(k) 202 1 fi exp [
(11)
Substituting Eq. (lob) into Eq. (11) and taking the ensemble probability density functidn over all values of k in the passband of the filter as the product of the individual probability functions yields the conditional probability density of S(k) given Am as
AWlAm)
=
’
(JSco)K
- &
Izd [S,(k)cos(2nkAm/M)
- s,(k)sin(2~kAmiM)12},
exp
(12)
where K is the number of spectral components contained within the passband. Maximizing this probability density function with respect to Am provides a maximum likelihood estimator for the shift. The resulting transcendental equation for the estimate, Am, is K-l C
k{2SR(k)S,(k)cos(4nkAh/M)
- [S;(k) - S:(k)]sin(4xkAfi/M)}*
= 0.
(13)
k=O
This equation can be solved by iterative techniques to determine the estimate for the shift. Once the shift has been determined, the spectrum of the data is multiplied by the phase term exp( - 27rrkAriz/M) to shift the estimated location of the axis of symmetry of the data back to the origin. Following this multiplication, the remaining imaginary components are by assumption necessarily due to noise, and so can be set equal to zero. This process completes the symmetrizing of the data, with the results that the inverse Hankel transform yields a real, even, and smoothed Abel inversion of the data. As a final step in the Abel inversion by means of transform methods, an error analysis based on the statistical properties of the input noise and the effective bandwidth of the applied filter can be made. It can be shown that for wide-sense-stationary input white noise with variance o* filtered with a bandwidth B and inverted via the Abel inversion integral in Eq. (2), the variance of the output noise is approximately’2 var[no(r)] N a2B/2r,
(14)
for sufficiently large r. If the imaginary portion of the noise is eliminated, as described previously, the output noise power is effectively halved and the factor in the denominator in Eq. (14) is 4 for the case of Abel inversion performed by the method presented here.
Abel inversion using transform techniques 4.
NUMERICAL
AND
371
EXPERIMENTAL
RESULTS
A number of numerical experiments have been carried out using both noise-free and noisy data to evaluate the performance of the integral transform method implemented using discrete transforms.‘O These numerical experiments using noise-free Gaussian and semi-elliptic data established the accuracy of the modified Candel method for the Hankel transform. Numerical experiments using noisy data were performed to establish the validity of the method used to symmetrize the data and to evaluate the effectiveness of the optimal filters used to smooth the data. A noisy trial function having a known inversion was generated by sampling a Gaussian function defined as Z,(X) = exp( -x*/400). The Abel inversion of Z,(x) is theoretically
(15)
given by’
6, (r) = (1/20&)exp(
- r*/400).
(16)
A zero mean sequence of random numbers having a Gaussian density function with a variance of 0.001 was generated and added to the sampled sequence of Eq. (15), which was then shifted by two integer sample points. Also, the sequence was zero padded so that the filtering by discrete Fourier transforms produces a linear convolution filtering. This trial sequence is shown in Fig. 1.
-.,ij -100
-
50
0
50
1
0
x
Fig. 1. Input noisy and off-center trial sequence I,(X).
-0.02
5 0
20
40
60
60
100
Fig. 2. Result of the Abel inversion for the sequence shown in Fig. I without filtering or symmetrizing. Circles mark the discrete calculated values connected by dotted lines while the theoretical inversion is shown by the solid line.
312
L. MONTGOMERY SMITHet al 0.04 r
I 20
-0.02 I 0
I 40
I 60
I 60
1 100
r Fig. 3. Result of the Abel inversion of the sequence shown in Fig. I with filtering and symmetrizing of the data. Circles mark the discrete calculated values connected by dotted lines, while the theoretical inversion is shown by the solid line.
The result obtained when this trial sequence was inverted without filtering or symmetrizing is shown in Fig. 2. The same trial function was then inverted using the symmetrizing procedure and an optimal filter having a passband from 0.0 to 0.0525 and a stopband from 0.0725 to 0.5. The result from this inversion is shown in Fig. 3. A principal advantage of this integral transform method is its computational efficiency. It has been applied to the inversions of a large number of images obtained from experiments with laser-sustained plasmas. I3 Each of these digitally sampled images consisted of 512 radial scans of the plasma and each radial scan consisted of 240 sample points. The inversion of 512 sets of 240-point data consumed about 8 min of CPU time when processed by the integral transform method, in contrast to a 20 h CPU time required for a curve-fit inversion using the method described in Ref. 6. The data analysis was carried out using a Masscomp MC500 Micro Supercomputer. The symmetrizing procedure generally consumed about 5-7 iterations for each set of data. Contour plots of the intensity and the emission coefficient for one of these images are given in Figs. 4 and 5, respectively. This particular image illustrates the performance of the method for intensity profiles which exhibit some asymmetry and which result in inversion profiles that are peaked on-axis (e.g. 49 mm), flat-topped (e.g. 45 mm) and having off-axis peaks (e.g. 43 mm).
3.0 2.4 1.8
E 2
1.2
0) 0.6 ,o i
2.4
t
Axial
Fig. 4. Contour
distance
(mm)
plot of a digitally sampled experimental intensity distribution laser-sustained plasma.
obtained
from a
373
Abel inversion using transform techniques
2.4 1.6 z 5
1.2
?j 0.6 ; .b! 0.0 v _ 0.6 0 G 1.2 B 1.6
3.0 40
43
46 Axial
distance
49
52
(mm)
Fig. 5. Contour plot of the plasma emission coefficient obtained by the Abel invesion of Fig. 4 using transform techniques.
5. CONCLUDING
REMARKS
A new integral transform method has been developed for the Abel inversion of measured intensity data in the presence of noise. The integral transform method has a number of advantages compared with existing curve-fit techniques. The curve-fitting methods do not allow for consideration of the spectral characteristics of the noise and the desired information or address the problem of symmetrizing the experimental data. The Abel inversion or the inverse Abel transform is obtained by performing a Fourier transform followed by the inverse Hankel transform. The efficiency of the well known FFT is exploited in the evaluation of both the Fourier and inverse Hankel transforms. This fast technique is made applicable to noisy and off-center experimental data by using frequency-domain filtering and a maximum likelihood estimator derived from the assumed noise characteristics. REFERENCES 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13.
H. R. Griem, Plasma Spectroscopy, McGraw-Hill, New York, NY (1964). 0. H. Nestor and H. N. Olsen, SIAM Rev. 2 2, 200 (1960). K. Bockasten, JOSA 51, 943 (1961). M. J. Freeman and D. Katz, JOSA 53, 1172 (1963). C. J. Cremers and R. C. Birkebak, Appl. Opt. 5, 1057 (1966). R. T. Shelby, Masters Thesis, Department of Mathematics, University of Tennessee, Knoxville, TN (1976). C. D. Maldonado, A. P. Caron, and H. N. Olsen, JOSA 55, 1247 (1965). R. N. Bracewell, The Fourier Transform and Its Applications, McGraw-Hill, New York, NY (1965). S. M. Candel, Comp. Phys. Commun. 23, 343 (1981). S. I. Sudharsanan, Masters Thesis, Department of Electrical Engineering, University of Tennessee, Knoxville, TN (1986). K. S. Shanmugam, Digital und Analog Communications Systems, Wiley, New York, NY (1979). L. M. Smith, IEEE Trans. Inform. Theory, to be published. R. D. Welle, D. R. Keefer, and C. Peters, AIAA Paper 86-1077, “Energy Conversion Efficiency in High-Flow Laser-Sustained Argon Plasmas,” Atlanta (1986).