This work builds on the theory laid out in [StefanovP2018] on sampling Fourier Integral Operators (FIOs). We discuss the specific application of Thermoacoustic Tomography, in which case the measurement operator is an FIO under suitable conditions. We discuss the theoretical resolution of given the sampling rate of and then discuss aliasing and averaged data. Lastly we will show empirical evidence of our findings using numerical simulations.
Thermoacoustic Tomography is a medical imaging method in which a short pulse of electromagnetic radiation is used to excite cells in some object we wish to image, typically the organs of a patient. Upon absorbing the EM radiation, the cells in the patient in turn vibrate, creating ultrasonic waves that then propagate out of the patient and are measured by any number of methods. Using this measured data, we then try to reconstruct, in some sense, an image of the inside of the patient. This is a hybrid imaging method which uses high contrast, low resolution EM radiation to excite the cells; and low contrast, high resolution ultrasound waves as measurement [Oraevsky1994, Kuchment2014, Kruger1999, Kruger2000, wang2015]. The hope is to be able to get an image with good contrast and resolution by combining these two types of waves.
More precisely, let be an open subset of Euclidean -space such that for some where is the Euclidean ball of radius . Suppose is a smooth function on supported in . We view as the initial pressure distribution internal to some object to be imaged. Then, after exposing to EM radiation, the ultrasonic waves created solve the acoustic wave equation:
Here, is the wave speed, which we take to be identically outside of . We assume that is a smooth function of . In addition, is the Riemannian metric on the space , assumed to be Euclidean on . We define , which is the metric form which determines the geometry of this problem. Assume is a solution to (1) for all . Further suppose that we have access to for where and is a relatively open subset of (for this paper, we will take ). We define for the distribution as the measurement operator:
where is the space of smooth functions on such that near . The methods used to collect data on are varied and include point detectors [Kunyansky2008, StefanovP2009, hristova2008, hristova2009], integrating line detectors [Burgholzer2006, Grun2007], circular integrating detectors [Haltmeier2007, Zangerl2009], and 2D planar detectors [Stefanov2017, Haltmeier2004]. We note that at least when
, by energy estimates,is well defined. We may actually even take to be a distribution in such that , and by conservation of energy, extends to a well defined operator. The closure of under the previously stated norm is the space , and we will assume unless otherwise stated.
1.1. as an FIO
To obtain an oscillatory integral representation of , we may use the geometric optics construction to solve for in up to a smooth error (see [StefanovP2009, Taylor81] for more details). This construction leads to the representation
where are solutions to the eikonal equation with initial conditions . Note that solutions to the eikonal equation are local in nature, and so this representation of is only valid until some time . However, we may then solve (1) with “initial” conditions and using the same geometric optics construction. In this way, we can obtain an “approximate” solution to (1) for all . Note by approximate, we mean up to a smooth error term. This error term could be quite large in the sense, but because it is a smooth term, it is negligible in the calculus of FIOs. It can be shown that is a sum of elliptic FIOs of order associated with locally diffeomorphic canonical relations that are each (locally) one-to-one mappings (see i.e. [StefanovP2018, StefanovP2009]). We record the canonical relations and here for later use:
Here, we have is the exit time of the geodesic starting at in the direction , is the point on the geodesic issued from at time and is the orthogonal (in the metric) projection of onto (the tangent bundle of the boundary of , so implicitly, we assume that is a at least a manifold). We assume that the metric induced by is non trapping, so that for all . Note that because each of the canonical relations and are one-to-one, the full canonical relation of the FIO given by is one-to-two, which makes intuitive sense as singularities split and travel along geodesics according to propagation of singularities theory.
The author would like to thank Dr. Plamen Stefanov for suggesting this problem and for his guidance in the analysis of this problem.
2. Preliminary definitions and theorems
2.1. Semiclassical analysis
The main definitions and theorems of semiclassical analysis and sampling that we use come from [Zworski2012, StefanovP2018]. For a more complete background on semiclassical analysis, see [Zworski2012]. In sampling the measurement operator
, we are interested in how the sampling rates affect our ability to resolve singularities with high frequency. To model this, we will rescale co-vectorsby a factor of where is a small parameter. We then examine families of functions (or distributions) that satisfy certain growth conditions as becomes small. Because of this, instead of considering the classical wave front set of a distribution, we consider the semiclassical wave front set, denoted . Note that is understood here to be a family of functions depending on the parameter , but we will drop this subscript when it will not cause confusion. A key tool in analyzing the behavior of the measurement operator
will be the semiclassical Fourier Transform, defined below.
Definition 2.1 (Semiclassical Fourier Transform).
The semiclassical Fourier transform of an -dependent family of distributions is defined as
If we denote the classical Fourier Transform by , then we have
Much like in classical analysis, we can use the semiclassical Fourier transform to define Sobolev norms on certain classes of functions or distributions.
Definition 2.2 (-Tempered family of distributions).
The -dependent family of distributions in is said to be -tempered if
is such that for some and . Here, we have .
Another key tool we will use is the idea of the semiclassical wave front set of an -dependent family of distributions.
Definition 2.3 (Semiclassical Wave Front Set).
The semiclassical wave front set of the -tempered family is defined to be the complement of the set of such that there exists with so that
for in a neighborhood of .
This set plays a similar role as the classical wave front set from microlocal analysis, however in general there is no sort of inclusion between these two sets. As an example [Zworski2012], the coherent state
has an empty wave front set in the classical sense, as it is a smooth function in both and , however its semiclassical wave front set is . Note also that the zero section is allowed to be a part of the semiclassical wave front set, unlike in the classical case. Also, we do not require the semiclassical wave front set to be a conic set, which is another way that this set differs from the classical wave front set.
We call elements of singularities, even though a function with finite semiclassical wave front set is actually smooth.
Definition 2.4 (-Do).
We will use the standard quantization to define semiclassical pseudodifferential operators. Fix and and let satisfy the following: For every and multi-indices and every compact set there exists some such that
for all and . We then say is a semiclassical symbol of order . Then we define the semiclassical pseudodifferential operator by
The -tempered family is said to be localized in phase space if there exists some such that
Note that because the functions we work with are semiclassically band limited (see definition 2.7), that all functions we work with can be assumed to be localized in phase space unless otherwise stated.
Definition 2.6 (Semiclassical Frequency Set).
For each tempered -dependent distribution localized in phase space, set
This is simply the projection of onto the second variable.
Definition 2.7 (Semiclassically Band Limited Functions).
We say that is semiclassically band limited (in ) if
is contained in an -independent set,
there exists a compact set such that for every open , we have for every there exists such that
Semiclassically band limited functions are those functions that can be reconstructed up to a smooth error from their samples, much like the band limited functions are those that can be perfectly reconstructed from their samples in the classical Nyquist Sampling theorem given a small enough sampling rate[Marks1991].
The main theorem used in [StefanovP2018] is the following:
Assume that , are open and bounded. Let satisfy
for some such that . Let be such that and near .
Assume that is an invertible matrix so that the images of
is an invertible matrix so that the images ofunder the translations , are mutually disjoint. Then for every ,
for every , and
The proof of this theorem essentially follows from the classical Nyquist sampling theorem and can be found in [StefanovP2018, Petersen1962]. For all applications in this paper, we take the matrix
above to be the identity matrix.
We make heavy use of the following theorem which relates how classical FIOs effect semiclassical wavefront sets from [StefanovP2018], where the reader can find the proof.
Let be an FIO in the class where is a Lagrangian manifold and . Then for every localized in phase space,
where is the canonical relation of .
This theorem shows how classical FIOs affect the semiclassical wavefront set away from the zero section. In particular, the semiclassical wavefront set of away from the zero section transforms in the same way the classical wavefront set does: it is transformed by the canonical relation associated with . The main assertion in [StefanovP2018] is that the sampling requirements of given are determined by , the canonical relation associated with .
3. Resolution limit of given sampling rate of
Suppose we wish to sample the at some fixed sampling rates and . Here we don’t assume that we know any information about , we only wish to see how fixing a sampling rate on affects our ability to resolve singularities of . Avoiding aliasing of is equivalent to (by Theorem 2.8)
where is the dual variable to , and is the dual variable to , with the th component of . Note that the norms and are taken in the corresponding metric. In particular, although the norm on is assumed to be Euclidean, the induced norm on the tangent space to the boundary, which we’ll call , is not necessarily Euclidean. We may use the canonical relation (2) associated with to write the inequalities above as
From this we can see that we have that avoiding aliasing is equivalent to
For most of the paper, we will assume that is Euclidean, although more general results hold.
3.1. The effect of on resolution
Consider the first inequality in (7) and assume that is taken small enough so as to not effect resolution of singularities of . The first inequality indicates that the sampling rate imposes a limit on the resolution of such that for fixed , there will be higher resolution of singularities of at points where the wave speed is slower, and likewise the resolution will be worse at those points where the wave speed is faster. In particular, given the relative sampling rate , we cannot resolve singularities at with frequency greater than
This is a local result. A global estimate for the maximum frequency of a singularity that is guaranteed to be resolved anywhere given the sampling rate is given by
3.2. The effect of on resolution
Assume now that is chosen small enough so as to not effect resolution of singularities of . The second inequality in (7)
tells us that the sampling rate imposes a limit on the resolution of such that singularities that intersect the boundary nearly perpendicularly will have higher resolution than those that hit the boundary nearly tangentially (at a large angle to the normal vector to at the point of intersection). Also, because is constant along the geodesic , we know in particular that where is the angle (in the metric) between and . This tells us that to avoid aliasing, we must have
We recall that , and in the case that is Euclidean, we get
For a fixed relative sampling rate , we cannot resolve singularities of of frequency greater than
Note in particular that if (i.e. the geodesic hits the boundary perpendicularly), then , and we will always be able to resolve the singularity at . Also note that this is a local result, and as is the case for “slow spots” in the speed give better resolution of singularities in general. Because , we also get the following estimate for the maximum frequency of a resolvable singularity, regardless of location:
Finally, because , we know , and we have the following (worst case) global estimate for the maximum frequency of a singularity of that can be resolved:
where is defined as before. In particular, we recover the result from [StefanovP2018] that for a semiclassically band limited with essential maximum frequency in the Euclidean case that we need to take sampling rates of satisfying
to avoid aliasing. These effects are shown in Figure 3.
3.3. CFL condition
We can relate this analysis to numerical solvers of the wave equation. When solving the wave equation numerically, a typical approach is to discretize the space and time domain, and use a finite difference scheme. Suppose we wish to simulate an experiment using a rectangular grid in the space coordinates and we collect data on the boundary of a square. Further, we assume that is Euclidean, and because the boundary is a rectangle, also the metric induced on the boundary is Euclidean. Suppose we have fixed each with a common value , where is the essential band limit on , i.e. . Note that by our choice of , there will not be aliasing of , provided is chosen well, as on the boundary in this rectangular grid, we have , where all of the as above have a common fixed step size . In order to choose , we recall that the frequency set is contained in the set . Because has a semiclassical band limit of , we know that , where is the projection onto the second factor. We know this because each . Also, by the analysis above, we know that , but . We also know that , so that the largest possible size of given the band limit on , is . It is then clear that we need to avoid aliasing. This tells us that we should take . Now, the CFL condition for the leapfrog finite difference scheme ([CFL59, Bartels2016, StrangNotes]) tells us that given a step size and wave speed , that we should take the time step to ensure stability of the finite difference scheme. But , because . This means, that if we’ve chosen , and we choose satisfying the CFL condition for the leapfrog finite difference scheme, then there will be no aliasing in the measured data at the boundary. Also, if , then the CFL condition is identical to the conditions on and required to avoid aliasing of the measured data .
4. Aliasing and artifacts
Now suppose that we know that is a semiclassically band limited function with essential band limit . In [StefanovP2018], it is shown that in order to avoid aliasing of , for a semiclassically band limited , we must have relative sample rates of and where is half the side length of a box bounding , is the sharp lower bound of the metric form on the unit sphere for all , and is the sharp upper bound on the induced metric on the Euclidean sphere in a fixed chart for . In the numerical examples that follow, is piecewise flat and parameterized in a Euclidean way, so that away from corners. Note that if is Euclidean, then setting , we have , and so that the relative sampling rates needed to avoid aliasing are
4.1. Under sampling in
. Then, by [StefanovP2018] there will be aliasing of . The error in the reconstruction can be modeled by the frequency shift operator
This operator is valid as long as (see Figure 4 (right)).
If we have not under sampled too critically in the variable, we would expect to only see this added error for , with more terms added as the under sampling becomes worse. As explained in [StefanovP2018], by Egorov’s Theorem, we expect to see artifacts in a reconstruction of that can be calculated by the canonical relation
where and can be calculated by finding the operator on the left. We do that now for :
where is the point of intersection of the geodesic issued from with , and where and . Aliasing artifacts are found using this mapping in Figures 5 and 6 below. The mapping is calculated in almost an identical fashion, however we have a change in sign in the variable.
We include a more complicated image reconstruction in Figure 16 along with the collected data in Figure 17. We also show how a smooth approximation of an line segment is affected by these artifacts in the image given in Figure 9. For this image and reconstruction, we have included the collected data and Fourier transform images in Figure 10.
4.2. Under sampling in
Now suppose that we have under sampled the variable, i.e. we have chosen for some . Then again, we will have aliasing and the error in the reconstruction will involve the frequency shift operator, but now will act on as
This operator is valid as long as . The canonical relation of the -FIO that operates on as a reconstruction of will then be given by (again, we only consider here)
where is the unit vector in the direction. Note that, in particular, this implies that the artifacts will have the same frequency as that of the original image, but perhaps with a space shift. Also, because this operator is valid as long as , if the geodesic emanating from hits the boundary perpendicularly, then the point will be unaffected by this shift in the reconstruction, i.e. there will be no artifacts that come from . This is true because if the geodesic emanating from hits perpendicularly, then and for any . Finding these artifacts in practice follows in much the same way as finding where artifacts occur for under sampling in the time variable. We illustrate this for the constant speed, Euclidean case in Figure 7 and see Figure 8 for the variable speed case.
We again include a more complicated image reconstruction in Figure 18 along with the collected data in Figure 19. We also show how a smooth approximation of an line segment is affected by these artifacts in the image given in Figure 11. For this image and reconstruction, we have included the collected data and Fourier transform images in Figure 12.
5. Averaged data
Suppose that the collected data has been averaged in the or variables for some reason (in practice this can be done to try to avoid aliasing, or in an attempt to reduce the noise in data). This can be modeled in a few ways, including taking a convolution with a smooth function that decreases away from the origin to 0. To model localized averaging however, we will consider data of the form , where is an -DO with a principal symbol of the form where is decreasing. The effect of is to limit , which will in principle remove the high frequency singularities of which will have a smoothing effect. From [StefanovP2018], we know that because is a FIO associated with the canonical map , that the composition can be written
where is a -DO with principal symbol where is the principal symbol of . So, for , we may calculate
Suppose we only average the time data in . This corresponds to taking above to give . This symbol takes its minimum values where
is maximized. Assuming for a moment thatis Euclidean, this means that we expect more blurring at points where the wave speed is “fast”. Additionally, we expect singularities with large frequencies to be blurred more than smaller frequencies where the wave speed is the same. These effects can both be seen in Figure 13.
Suppose now that we only average data in the spatial variable . This corresponds to taking above and we get the principle symbol of to be
Here the norm is the induced norm on the boundary, which we have noted in this paper as . This symbol takes its smallest values when is large, i.e. when the geodesic issued from intersects the boundary at a large angle. In addition, we expect singularities that hit the boundary perpendicularly to be affected far less by averaging of data in the variable. In addition, because where is the angle between and we expect to see more blurring at points with faster speeds or higher frequency. For constant speeds , the effect of averaging data in is uniform in , but the effect is local for averaging in , due to the blurring depending on the angle of intersection made by geodesics. In addition, with a variable speed singularities in “slow spots” of will have higher resolution when blurring in the -data, but their resolution will still depend on how geodesics hit the boundary. The result is a roughly uniform blurring in fast spots of