Sampling in Thermoacoustic Tomography

12/09/2019 ∙ by Chase Mathison, et al. ∙ 0

We explore the effect of sampling rates when measuring data given by Mf for special operators M arising in Thermoacoustic Tomography. We start with sampling requirements on Mf given f satisfying certain conditions. After this we discuss the resolution limit on f posed by the sampling rate of Mf without assuming any conditions on these sampling rates. Next we discuss aliasing artifacts when Mf is known to be under sampled in one or more of its variables. Finally, we discuss averaging of measurement data and resulting aliasing and artifacts, along with a scheme for anti-aliasing.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 6

page 10

page 12

page 13

page 14

page 17

page 19

page 20

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. Introduction

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:

(1)

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:

(2)

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.

Acknowledgments

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-vectors

by 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

Definition 2.5.

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

  1. is contained in an -independent set,

  2. is tempered,

  3. 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].

2.2. Sampling

The main theorem used in [StefanovP2018] is the following:

Theorem 2.8.

Assume that , are open and bounded. Let satisfy

(3)

for some such that . Let be such that and near .

Assume that

is an invertible matrix so that the images of

under the translations , are mutually disjoint. Then for every ,

(4)

for every , and

(5)

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.

Theorem 2.9.

Let be an FIO in the class where is a Lagrangian manifold and . Then for every localized in phase space,

(6)

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

(7)

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

(8)

This is illustrated in Figures 1 and 2 below.

Original Image
Time Reversal Reconstruction
Wave Speed with Fast Spot
-2
2
-2
2
-2
2
-2
2
-2
2
-2
2
Figure 1. Resolution of given a fixed sampling rate of . The wave speed here has a fast spot centered at . We can see that this is precisely where the reconstruction of has poor resolution when under sampled in the variable, as explained above.
Original Image
Time Reversal Reconstruction
Wave Speed with Slow Spot
-2
2
-2
2
-2
2
-2
2
-2
2
-2
2
Figure 2. Resolution of given a fixed sampling rate of . The wave speed here has a slow spot centered at . We can see that this is precisely where the reconstruction of has the best resolution when under sampled in the variable, as explained above.

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:

(9)

We note that if one wants to be able to resolve singularities of with frequency , then by considering (8) and (9), the sampling rates and of should be taken to be at least

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 .

Original Image
Time Reversal Reconstruction
Wave Speed with Fast Spot
-2
2
-2
2
-2
2
-2
2
-2
2
-2
2
Figure 3. Resolution of given a fixed sampling rate of the space variables on the boundary . We can see that the blurring effect is roughly uniform for points near the fast spot in the wave speed , but that there are singularities in the region where far from the fast spot that are also highly affected. These singularities hit the boundary with a larger angle to the outward pointing normal vector, and so we expect lower resolution there.

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)).

-
Figure 4. The characteristic cone in which must lie. The cone on the left shows the possible range of the covector which is determined by the canonical relation associated with . The image on the right shows the possible range of covectors after under sampling (in ). Note that the red regions have been shifted up and down from the original frequency set by translation due to under sampling.

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.

Original Image
Reconstructed Image
Figure 5. Tracing the aliasing artifacts by using geodesics. We have used the constant wave speed for this example. Here we have under sampled in and show the image of the singularity under the canonical relations given by for . Note that the low frequency singularity does not cause artifact, but the high frequency singularity vanishes in the reconstruction and causes aliasing artifacts.
Original Image
Reconstructed Image
Figure 6. Artifacts in a reconstructed image with under sampled in time variable and a variable wave speed. We trace the geodesics to find the image of under the map as explained above.

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.

Original Image
Reconstructed Image
Figure 7. Artifacts in a reconstructed image with under sampled in space variables. Here we take . Specifically, here was under sampled on the left and right edges of the square. Note that there is no artifact in the reconstructed image coming from the pattern in the upper right corner of the square, because singularities from this pattern hit the boundary of the square perpendicularly. Note also that the original singularity still remains with half its amplitude because we did not under sample along the bottom edge of the square.
Original Image
Reconstructed Image
Figure 8. Artifacts in a reconstructed image with under sampled in space variables and a variable wave speed. Specifically, here was under sampled on the top and bottom edges of the square. The artifacts in the reconstruction have the same frequency as the original, but with a space shift due to under sampling.

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.

Original Image
Reconstructed Image
Figure 9. Original and reconstructed image of a smooth approximation of an line segment. Here we have under sampled in . The under sampling has resulted in blurring of this “line segment”. This is due to the fact that under sampling in shifts high frequency data in .
Collected Data
Downsampled Data
Fourier Transform
of Collected Data
Fourier Transform of
Undersampled Data
Figure 10. Collected data and Fourier transform along with under sampled data in for example given in Figure 9. Data was collected on all edges of the square at a rate guaranteeing no aliasing. Shown is the data from the bottom edge of the square. We can see that under sampling in has resulted in the Fourier Transform of being folded into the band limit region. Under sampling in shifts large frequencies from , thus producing the blurred image we see in the right of Figure 9.
Original Image
Reconstructed Image
Figure 11. Original and reconstructed image of a smooth approximation of an line segment. Here we have under sampled in . This has resulted in some blurring, but also in high frequency artifacts.
Collected Data
Downsampled Data
Fourier Transform
of Collected Data
Fourier Transform of
Undersampled Data
Figure 12. Collected data and Fourier transform along with under sampled data in for example given in Figure 11. In contrast to when we under sample in , we see that high frequencies in are not necessarily eliminated when we under sample in , but there is a phase shift. This results in more high frequency artifacts in the image on the right in Figure 11.

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 that

is 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.

Original Image
Wave speed with fast spot
Reconstructed image with
data averaged in variable.
Figure 13. Reconstructed image from data that has been averaged in time variable. We can see that the reconstructed image is most blurred at the points where the speed is fast, and there is less blurring where .

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