The derivative-free Fourier shell identity for photoacoustics

In X-ray tomography, the Fourier slice theorem provides a relationship between the Fourier components of the object being imaged and the measured projection data. The Fourier slice theorem is the basis for X-ray Fourier-based tomographic inversion techniques. A similar relationship, referred to as the ‘Fourier shell identity’ has been previously derived for photoacoustic applications. However, this identity relates the pressure wavefield data function and its normal derivative measured on an arbitrary enclosing aperture to the three-dimensional Fourier transform of the enclosed object evaluated on a sphere. Since the normal derivative of pressure is not normally measured, the applicability of the formulation is limited in this form. In this paper, alternative derivations of the Fourier shell identity in 1D, 2D polar and 3D spherical polar coordinates are presented. The presented formulations do not require the normal derivative of pressure, thereby lending the formulas directly adaptable for Fourier based absorber reconstructions.

component of the pressure data and its normal derivative to a specified spatial Fourier component of the source object function. The difficulty with Anastasio et al's derivation is that knowledge of the normal derivative of pressure is required-a value that is not normally obtained through experiments.
In this paper, the Fourier shell identity in 1, 2 and 3 dimensions, such that the normal derivative is not required, are derived and presented. This provides a direct relationship between the temporal Fourier components of the pressure data and the Fourier components of the object function. The Fourier shell identity as derived in this paper requires no additional information than that typically acquired in experiments and thus can be directly used for absorber reconstructions.

Photoacoustic governing equations
Diebold gives a concise explanation of the governing equation for the pressure that results from launching a photoacoustic wave (Diebold 2009). The governing equation is given by where β is the thermal expansion coefficient, c s is the speed of sound, C p is the specific heat, H is the energy per unit volume and time deposited by the optical radiation beam, and p( r, t) is the pressure of the acoustic wave, a function of space ( r) and time, t. As is common, it is assumed that H is a separable function of space and time, so that H (� r, t) = A(� r)I(t). In this work, the temporal Fourier transform is employed with the angular frequency, non-unitary convention for forward and inverse transforms, as given by where a tilde (∼) over the variable has been used to denote a temporal Fourier transform, and assuming suitability of the function f for Fourier transformation. The same formulation can be used to define spatial Fourier transforms, where an overhat notation is used to denote a spatial Fourier transform where the spatial variable r is transformed to the spatial frequency variable ω, so for example f ( r) is spatially Fourier transformed to f ( � ω). Taking the temporal Fourier transform of (1), it then follows that where k = ω/c s is the wave number. Anastasio et al. (2007) demonstrated a mathematical relationship between the pressure wavefield data function and its normal derivative measured on an arbitrary aperture that encloses the object and the three-dimensional (3D) Fourier transform of the optical

Fourier shell identity
absorption distribution evaluated on concentric spheres. They referred this relationship as a 'Fourier-shell identity' , analogous to the well-known Fourier slice theorem of X-ray tomography. The Fourier-shell identity as derived by Anastasio et al. (2007) is given by where Ω 0 is an arbitrary measurement surface that is smooth, closed and encloses the object, dS ′ is the differential surface element on Ω 0 , n ′ is the outward normal vector to Ω 0 at the point � r ′ 0 ∈ Ω 0 , ŝ is an arbitrary unit vector and the Â ( � ω) indicates the 3D spatial Fourier transform of A( r).
In practical applications, the normal derivative of the pressure wavefield will not be measured so the Fourier shell identity as given in (4) is not immediately useful. Equation (4) permits determination of concentric 'shells' of Fourier components of Â � ω = kŝ from knowledge of p � r ′ 0 , k and its derivative along the n ′ . Because ŝ can be chosen to specify any direction, Â � ω = kŝ specifies the Fourier components of Â ( � ω) that reside on a spherical surface of radius |k|, whose center is at the origin.
It is shown in the proceeding sections how Eq. (4) can be made specific for the 1D, 2D polar and 3D spherical polar cases. The identity is also re-derived in those three cases so that the derivative of pressure does not appear in the formulation.

Fourier shell identity in 1D
Using the same approach as used for the original derivation in Anastasio et al. (2007), it is shown in the "Appendix" that the 1D expression of the Fourier shell identity is given by Equation (5) is the 1D equivalent of Eq. (4), the Fourier shell identity. In Eq. (5), the overhat indicates a spatial Fourier transform in the 1D spatial variable, z, meaning that the spatial variable z transforms to the spatial frequency variable ω z . Equation (5) at first appears different from Anastasio's formulation in Eq. (4). The source of this apparent difference is that in it is necessary to enclose the source function and to do so in 1D, the "detecting surface" must be located at both z = ±z 0 . This follows because a detector at only one of z = ±z 0 will not enclose the source.

Multivariable Fourier analysis
To proceed with the derivation in this manuscript, the photoacoustic Eq. (1) is considered in one dimensional Cartesian coordinates so that position is a function of z only, � r = (z). Taking the temporal Fourier transform of (1) (denoted with an over-tilde) so that time, t, transforms to the temporal frequency variable ω, and then further taking a spatial Fourier transform (denoted with an overhat) in the spatial variable, z, so that spatial variable z transforms to the spatial frequency variable ω z , leads to where k = ω/c s is the wave number.

The Fourier shell theorem in 1D
From (6), the inverse spatial Fourier transform can be computed with the help of the appropriate choice of Theorem 5 from Baddour (2011). The third option of Theorem 5 from Baddour (2011) states that the following result holds true: Assuming that Â (ω z ) has no poles and remains bounded, then the application of the inverse spatial Fourier transform (2) to (6), while making use of the identity in (7) gives Equation (8) is the statement of the Fourier shell identity in 1D. This result is very powerful as it relates the temporal Fourier components of pressure p(z, ω) directly to Â (k), which represents the Fourier components of the spatial Fourier transform of the absorber evaluated at ω z = k where k = ω/c s . It is shown in the next section that this is exactly the same as the result derived by Anastasio et al. (2007) for the Fourier shell identity, however Eq. (8) does not require knowledge of the derivative of pressure and thus can be easily inverted. For example, Eq. (8) implies measurements made at some detector location z = z 0 < 0 (measurement in reflection) can be used to write Hence, the success in being able to reconstruct A(z) is partly in the ability to gather sufficient (k intensive) information about Â (k) to enable the inverse spatial Fourier transform to be computed.
Assuming z = z 0 < 0 (which corresponds to a measurement in reflection) and using (9), then Equation (10) requires measurements of p(z 0 , ω) and Ĩ (ω) at a sufficiently wide bandwidth to enable the Fourier integral to be computed or approximated. In Eq. (10), the z 0 /c s term in the exponential represents the time taken for the signal to travel from the absorber to the detector and that time is an indication of the spatial location of the absorber. The term in front of the integral is simply a scaling term. Therefore, Eq. (10) states that the absorber profile in space can be found by finding the inverse temporal Fourier transform Ā (t) = F −1 p(z 0 , ω)/Ĩ(ω)|ω → t via well known numerical techniques for calculating Fourier transforms and then scaling A(z) =Ā(c s t) to find the shape of the absorber function as a function of space.

Comparison with the Fourier shell identity
Now, Eq. (5) is the "Fourier shell identity" as put forward by Anastasio et al. (2007), where the temporal Fourier component of pressure values p(z 0 , k) and its derivative are related to the values of a specified Fourier component, Â (k). The problem with the formulation in (5) is the requirement for measurements of the derivative of p(z 0 , k), which does not occur in practice. However, it can be shown that the formulation presented in this paper, as given in Eq. (8), leads to the same result in Eq. (5).
Using the newly derived form of the Fourier shell theorem as given in Eq. (8), evaluate Simplifying (11) gives Hence, the version of the Fourier shell theorem in 1D as presented in Eq. (8) satisfies the Fourier shell identity in 1D as derived by Anastasio et al. (2007) and shown in Eq. (5). However, the formulation of Eq. (8) provides a useful alternative since no knowledge or measurement of the derivative of pressure is required.

2D polar coordinates
For a 2D function in polar coordinates, the function can be written as f (� r) = f (r, θ) and the θ dependence can be expanded into a Fourier series due to the 2π periodicity of the function in θ. Hence, where the Fourier coefficients f n (r) can be found from The process of finding f n (r) from f (r, θ) can be thought of as a "forward Fourier series transform" where the continuous θ variable is transformed to the discrete (although infinite) n variable. The reverse transform is the process of performing the summation over the discrete n variable, as given in Eq. (13), to recover the original function as a continuous function of θ. The use of the Discrete Fourier Transform (DFT) via fast algorithms such as the FFT (Cooley and Tukey 1965) to numerically compute the continuous, infinite expressions of (13) and (14) has received much attention in the literature. There is a large body of work on how well the continuous Fourier transform or series can be approximated with the discrete FFT, for example in Epstein (2005). A transform that is quite useful in 2D polar coordinates is the Hankel transform of order n, defined by the integral (Chirikjian and Kyatkin 2000) where J n (z) is the nth order Bessel function. In 2D coordinates, the overhat symbol f is used to denote a Hankel transform. The superscript (n) is used as a reminder of the order of the Hankel transform-in this case an nth order Hankel transform. If n > −1/2, the transform is self-reciprocating and the inversion formula is given by The Laplacian in 2D polar coordinates is given by For a 2D function in polar coordinates as given by (13), the Laplacian gives Hence, for a function written in the form of Eq. (13), the required form of the Laplacian is denoted with ∇ 2 n where this operator is defined by Consider the application of an nth order Hankel transform to the 2D Laplacian ∇ 2 n of any function of r only, H n ∇ 2 n g(r) : where the second line follows from a simple application of integration by parts and the definition of the nth order Hankel transform. In other words, the nth order Hankel transform of the nth order Laplacian of any function g(r) is the product of −ρ 2 and the nth order Hankel transform of g(r). As is familiar from Fourier theory, derivatives transform to products under a suitable choice of integral transform. It is noted that the use of the Laplacian ∇ 2 n necessitates an nth order Hankel transform to give the desired simple result in (20).

Fourier shell identity in 2D
To proceed with the derivation in 2D, the equation for pressure (3) is again considered, and only the forward Fourier series transform is taken. Specifically, p n (r, k) are the Fourier series coefficients of p(r, θ , k) via (14) and similarly A n (r) are the Fourier series coefficients of A(r), also via Eq. (14). Hence, from (3), the relationship between them is given by It is shown in the "Appendix" that the Fourier shell identity stated in 2D is given by In Eq. (22), Â (n) n (k) is the nth order Hankel transform of A n (r) evaluated at ρ = k and r = r 0 is the location where a measurement is made. The tilde refers to a temporal Fourier transform, the overhat refers to a forward Hankel transform and the superscript (n) indicates the order of the transform. Equation (22) is the 2D statement of Eq. (4), the Fourier shell identity, where it is noted that the value of p n (r 0 , k) and its derivative in the radial direction are required.

Multivariable Fourier analysis
From the equation for pressure (3), two dimensional polar coordinates are assumed so that position is a function of radius and angle, � r = (r, θ ). Then, the forward Fourier series transform it taken, followed by an nth order Hankel transform. Specifically, p(r, θ , k) transforms to p n (r, k) via (14), which in turn transforms to p (n) n (ρ, k) via (15), an nth order Hankel transform of the nth term in the Fourier series for p(r, θ , k). Symbolically, p(r, θ, k) →p n (r, k) →p (n) n (ρ, k). The same applies to the absorber function so that A(r, θ ) → A n (r) →Â (n) n (ρ), and furthermore the Laplacian transforms under the same set of operations as ∇ 2 → ∇ 2 n → −ρ 2 . The overhat with superscript (n) denotes an nth order Hankel transform, and the subscript of n indicates the nth term of the Fourier series. Taking transforms and then cleaning up gives

Fourier shell theorem in 2D polar coordinates
It is shown in Baddour (2011) that the following result holds true Here, φ is an analytic function defined on the positive real line that remains bounded as x goes to infinity (has no poles), H (2) n (x) is a Hankel function of order n. Given the definition of the Fourier transform that is being currently used, the result in (24) satisfies the Sommerfeld radiation condition, ensuring an outwardly propagating wave.
Returning to Eq. (23), the inverse Hankel transform can be taken the help of the theorem shown in Eq. (24). Assuming that Â (n) n (ρ) has no poles and remains bounded, then Equation (25) is the statement of the Fourier shell theorem in 2D. This result if very powerful as it relates the temporal Fourier components of pressure p n (r, k) to the Hankel components of Â (n) n (k). Equation (25) is also easily invertible. Meaning, given measurements made at some detector location r = r 0 , then Hence, the success in being able to reconstruct A(� r) = A(r, θ ) is partly in the ability to gather sufficient (k intensive) information about Â (n) n (k) to enable the inverse Hankel transform to be computed via Equations (26) and (27) require measurements at a sufficiently wide bandwidth to enable the integral in (27) to be computed or approximated. The second challenge is the computation of a sufficient number of A n (r) terms (sufficient n) to enable A(� r) = A(r, θ ) to be reconstructed from (or an approximation to) From (27), a sufficient number of A n (r) terms is a consequence of a sufficient number of p n (r 0 , k) terms, which means that sufficient angular data must be collected to enable p n (r 0 , k) to be calculated or approximated from

Comparison with the Fourier shell identity
Equation (22) is the Fourier shell identity as put forward by Anastasio et al. (2007), where the temporal Fourier component of pressure values p n (r, k) and its derivative are related to the values of a specified Hankel component, Â (n) n (k). The problem with the formulation in (22) is the requirement for measurements of the derivative of p n (r, k), which does not occur in practice. However, it can be shown that the present formulation, as given in Eq. (25) leads to the same result as in (22). Using the expression for p n (r, k) given in (25), it can be substituted into the left hand side of (22) to evaluate which simplifies to There exist well known Wronskian relationships for the Bessel functions (Abramowitz and Stegun 1964;Olver et al. 2010), given by Hence, using Eq. (32), Eq. (31) simplifies to This is exactly the same expression as on the right hand side of (22), which implies that the results derived in the previous section and those derived by Anastasio et al. (2007) are identical. However, the formulation proposed here, namely Eq. (25) provides a useful alternative since no knowledge or measurement of the derivative of pressure is required.

3D spherical polar coordinates
For a 3D function in spherical polar coordinates, f (� r) = f (r, ψ, θ ), where 0 ≤ ψ ≤ π represents the colatitude which ranges from 0 at the North Pole, to π/2 at the Equator, to π at the South Pole and 0 ≤ θ ≤ 2π represents the longitude (azimuth angle). The angular dependence can be expanded into a spherical harmonic series (similar to a Fourier series in 2D polar coordinates) so that the function can be written as where the spherical harmonic coefficients in the series are given by The Y m l (ψ, θ) are spherical harmonics, solutions to the angular portion of Laplace's equation in spherical polar coordinates, and can be shown to be orthogonal. These spherical harmonics are given by (Chirikjian and Kyatkin 2000) where Y m l is a called a spherical harmonic function of degree l and order m, and P m l is an associated Legendre function. The coefficients f m l (r) are sometimes referred to as the spherical Fourier transform of f (r, ψ, θ ) ). The process of finding f m l (r) from f (r, ψ, θ ) can be thought of as a "forward spherical Fourier series transform" where the continuous (ψ, θ ) variables are transformed to the discrete (although infinite) m, l variables. The reverse transform is the process of performing the summation over the discrete m, l variables, as given in (34), to recover the original continuous function. A transform that is useful in spherical polar coordinates is the spherical Hankel transform and its inverse. These are defined as (Bracewell 1999;Piessens 2000) S n is used to specifically denote the spherical Hankel transform of order n. The spherical Hankel transform is particularly useful for problems involving spherical symmetry.
The Laplacian in 3D spherical polar coordinates is given by For a 3D function in spherical coordinates as given by (34), the Laplacian simplifies to Hence, from (39) it can be seen that for a function written in the form of (34), the required form of the 3D Laplacian is denoted with ∇ 2 l where this operator is defined by Consider the application of an nth order spherical Hankel transform to the 3D Laplacian ∇ 2 n of any function of r only, S n ∇ 2 n g(r) : (41) S n ∇ 2 n g(r) = ∞ 0 d 2 g(r) dr 2 + 2 r dg(r) dr − n(n + 1)g(r) r 2 j n (ρr) r 2 dr A simple application of integration by parts along with the definition of an nth order spherical Bessel function gives In other words, the nth order spherical Hankel transform of the nth order 3D Laplacian of g(r) is the product of −ρ 2 and the nth order spherical Hankel transform of g. As seen before, derivatives transform to products under the appropriate choice of transform.

Fourier shell theorem in 3D spherical polar coordinates
The analysis returns to the equation for pressure, (3) and proceeds by taking only the forward spherical harmonic transform. Specifically, p(r, ψ, θ, k) transforms to p m l (r, k) via (35), the Laplacian transform to ∇ 2 l and finally A(r, ψ, θ ) transform to A m l (r), also via Eq. (35). This gives It is shown in the "Appendix" that the Fourier shell identity stated in 3D is given by Equation (44) is 3D statement of Eq. (4), the Fourier shell identity. In (44), the tilde refers to a temporal Fourier transform, the overhat refers to a forward spherical Hankel transform, and the (l) superscript refers to order of the spherical Hankel transform.

Summary and conclusions
In this paper, re-derived and derivative-free formulations of the Fourier shell identity were presented in 1D, 2D polar and 3D spherical polar coordinates. These were shown to be equivalent to the previously derived Fourier shell identity (Anastasio et al. 2007), however knowledge of the derivative of pressure is not required, hence lending the formulas directly applicable to Fourier-based absorber reconstruction schemes. The formulas are restated here as a summary. In 1D, the Fourier shell identity is given by In Eq. (57), the ∼ refers to a temporal Fourier transform and the overhat refers to a 1D spatial Fourier transform.
In 2D polar coordinates, the Fourier shell identity is given by In Eq. (58), the n subscript refers to the nth Fourier coefficient in the Fourier series. The overhat refers to a forward Hankel transform, the superscript (n) refers to the order of the Hankel transform and H (2) n (x) is a Hankel function of order n. In 3D spherical polar coordinates, the Fourier shell identity is given by In Eq. (59), the m l subscript/superscript refer to the m l term in the spherical harmonic series expansion, the overhat refers to a forward spherical Hankel transform, the (l) superscript refers to order of the spherical Hankel transform and h (2) n (x) is a spherical Hankel function of order n.
It was assumed in Anastasio et al. (2007) that the measurement surface Ω 0 , which is here r = r 0 , completely encloses the source function. This means that A(r, ψ, θ ) has compact support and the integration limits on the right hand side of (82) can be taken to infinite without loss of generality. It then follows that where the overhat indicates a spherical Hankel transform and the Hankel transform