Semiclassical Sampling and Linear Inverse Problems Sampling meets - - PowerPoint PPT Presentation

Semiclassical Sampling and Linear Inverse Problems

Plamen Stefanov

Sampling meets Microlocal Analysis

IAS 2019

The classical Shannon sampling theorem says that if 𝑔𝑦 has a Fourier transform 𝑔 supported in the box 𝐶, 𝐶 (i.e., it is band‐limited), then 𝑔𝑦 is uniquely and stably determined by its samples 𝑔 𝑡𝑙 , 𝑙 ∈ 𝒂 if the sampling rate 𝑡 satisfies 0 𝑡 𝜌/𝐶. More precisely, 𝑔 𝑦 𝑔 𝑡𝑙 𝜓 𝜌 𝑡 𝑦 𝑡𝑙

• ∈𝒂

, 𝜓 𝑦 sinc𝑦

• and (unitarity)

𝑔 𝑡 𝑔 𝑡𝑙

• ∈𝒂

. The proof is simple. Think of 𝑔 as the inverse FT of 𝑔 . Then the samples 𝑔𝑡𝑙 are the Fourier coefficients of 𝑔 , more precisely of its periodic extension over the lattice 𝐶𝒂.

Graph of sinc The red curve is the function reconstructed, indistinguishable from the

• riginal. The blue curve is a cubic spline interpolation. Note that the

reconstruction has a global character. Here 𝑡 0.95 ∗ Nyquist.

In the frequency domain, things look like this: Nyquist condition satisfied (no aliasing). Since 𝐶, 𝐶 ⊂ 𝜌/𝑡, 𝜌/𝑡, 𝜌 𝑡 𝜌 𝑡 3𝜌 𝑡 3𝜌 𝑡 𝑔

• 𝟐/,/

𝜓

• Take the 2𝜌/𝑡‐periodic extension 𝑔
• f 𝑔

; then the samples are its Fourier coefficients: 𝑔

• 𝜊 𝑡 𝑔 𝑡𝑙 𝑓
• .

Multiply that extension by 𝟐 /,/ 𝜊 to get 𝑔 back: 𝑔 𝜊 𝟐 /,/ 𝜊 𝑡 𝑔 𝑡𝑙 𝑓.

• Take ℱ. The effect of the multiplication is the appearance of ℱof

𝟐 , 𝜊 in a convolution, which is the sinc function. If we have oversampling (supp𝑔 strictly in 𝐶, 𝐶), we can choose 𝜓̂ smooth, hence 𝜓 will be in the Schwartz class unlike the sinc function. 𝐶 𝐶

Aliasing: If the Nyquist condition is not satisfied: 𝐶 𝐶 3𝐶 3𝐶 𝑔

• 𝟐,

The 2𝐶‐periodic extension has overlapping segments. If we restrict back to 𝐶, 𝐶, we get 𝑔 truncated plus other stuff, i.e., we get 𝑔 modulo 2𝐶. Its inverse FT is not 𝑔 anymore. This creates aliasing. Frequencies get shifted. If we have a smooth 𝜓̂ instead of 𝟐,, we get an FIO, actually (in the sense below).

Oversampled

• n 81x81 grid

Undersampled

• n 41x41 grid

|𝑔 | The FT of the undersampled 𝑔 ℱ ℱ Aliasing in 𝐒: The original consist of two patterns: one higher frequency than the other. First, we sample and reconstruct properly. Next, we undersample the higher frequency pattern but still sample properly the lower frequency one. The reconstruction changes the direction and the frequency of the undersampled pattern. The Fourier transforms demonstrate the shifting (folding) of the frequencies.

We are interested in sampling the data g 𝐵𝑔 with 𝐵 linear, 𝑕 given, 𝑔 ? We assume that 𝐵 is an FIO, which is true very often: in integral geometry, thermo‐ acoustic tomography, etc. (i) Sampling 𝑩𝒈: Having an estimate of supp 𝑔 (the smallest detail of 𝑔), how dense should we sample 𝐵𝑔? (ii) Resolution limit on 𝒈, given the sampling rate of 𝑩𝒈. (iii) Aliasing: if we undersample 𝐵𝑔, what kind of aliasing we get for 𝑔? (iv) Averaged measurements/anti‐aliasing: if we blur the data (either because the detectors are not points and average already or because we want to avoid aliasing), what do we get for 𝑔? To answer (i), one may say – estimate supp 𝐵𝑔 knowing supp 𝑔 and apply the sampling theorem, we are done. The problem is that this is not

• straightforward. 𝐵𝑔 may not be even band limited if 𝑔 is.

We look at the problem as an asymptotic one. The “smallest detail” is a small parameter tending to 0. Rescale 𝜊 to 𝜊/ℎ (i.e., 𝜊 𝜃/ℎ) with 0 ℎ ≪ 1. This makes it a semi‐classical problem! If 𝜊 𝐶/ℎ, then for the rescaled 𝜊 we have 𝜊 𝐶. We want to work in the phase space, i.e., to account for the 𝑦 dependence as well. So band limited now means that 𝑔

𝑦 has a compact WF 𝑔 .

(v) Semi‐classical sampling. How to sample 𝑔

𝑦 with a compact WF 𝑔 ?

• Uniform sampling on rectangular or non‐rectangular lattices?
• Non‐uniform sampling?
• How many sampling points are enough?

This brings us to the first problem we need to study:

Theorem (semi‐classical sampling): Let 𝑔

∈ 𝐷 Ω with WF 𝑔 ⊂ Ω 𝐶, 𝐶 . Then for 𝑡 𝜌/𝐶,

𝑔

𝑦

𝑔 𝑡ℎ𝑙 𝜓 𝜌 𝑡ℎ 𝑦 𝑡ℎ𝑙

∈𝒂

𝑃𝒯 ℎ , where 𝜓 is a product of sinc functions. Parseval’s equality holds, too, up to 𝑃 ℎ . Just a rescaled classical version, with error estimates. The condition on WF 𝑔 is a condition on ℱ𝑔 modulo 𝑃 ℎ . The step size is 𝑡ℎ with 𝑡 𝜌/𝐶. As above, if WF 𝑔 ⊂ Ω 𝐶, 𝐶 (oversampling), 𝜓 can be made rapidly decreasing. We call the projection Σ𝑔 of WF 𝑔 onto 𝜊 the frequency set of 𝑔. (v) Semi‐classical sampling

𝐶 𝐶 𝐶 𝑡ℎ 𝜌ℎ 𝐶 𝑔 𝑔

• The sampling lattice

The Fourier domain 𝐶 Σ𝑔 (v) Semi‐classical sampling

𝐶/2 𝑔 The sampling lattice The Fourier domain The lattice could be rectangular 𝐶/2 𝐶 𝑔

• 𝐶

Σ𝑔 𝐶/2 (v) Semi‐classical sampling 𝑡ℎ 𝜌ℎ 𝐶

𝐶 𝐶 𝐶 𝑔

• The sampling lattice

The Fourier domain The previous sampling lattice (rescaled by 2) works for Σ𝑔 supported like

• this. The reason is that we can tile the plane with (non‐intersecting) shifts

generated by the vectors 𝐶/2,0 and 0, 𝐶. Then 𝜓̂ has to be supported there. ‐𝐶/2 𝐶/2 Σ𝑔 (v) Semi‐classical sampling 𝑔 𝑡ℎ 2𝜌ℎ 𝐶

One can use linear transformations. If shifts of Σ𝑔 under the translations 𝜊 ↦ 𝜊 2𝜌𝑋∗ do not intersect, where det 𝑋 0, then we can sample over the lattice 𝑦 𝑡ℎ𝑋𝑙, 𝑙 ∈ 𝒂 with 𝑡 1. This is a periodic but not a rectangular lattice (a parallelogram one). One can reduce the number of samples by taking some partition of unity 𝜓 and in each one, sampling according to Σ𝜓𝑔. This gets us closer to the important question of non‐uniform sampling. In the classical case, there is a well cited (780 google scholar citations) paper by LANDAU in Acta, saying, among the rest, that if we can sample uniquely and 𝑀‐ stably some 𝑔 over a possibly non‐uniform set of points then the density of that set must have a Weyl type of lower bound proportional to the Lebesgue measure of supp 𝑔 . This links sampling to spectral theory. (v) Semi‐classical sampling

We have the following microlocal analog of this result.

• Theorem. Let 𝑦ℎ, 𝑘 1, … 𝑂ℎ be a set of points in 𝑺. Let 𝐿 ⊂

𝑈∗𝑺 be a compact set. If 𝑔 𝐷 𝑂 ℎ 𝑔

𝑦 ℎ

• 𝑃 ℎ

for every 𝑔

with WF 𝑔 ⊂ 𝐿, then

𝑂 ℎ 2𝜌ℎ Vol𝐿1 𝑝 ℎ . The main novelty here, besides the semi‐classical setting, is that we relate the number of points needed to the phase volume of 𝑔

, i.e., to

Vol WF 𝑔 instead of Volsupp 𝑔 (multiplied by Volsupp 𝑔 ). The proof uses the spectral asymptotics (DIMASSI & SJÖSTRAND) for the smoothened version of the semiclassical ΨDO 𝟐𝑦, 𝜊. In the uniform cases above, we have equality. (v) Semi‐classical sampling

(i) Sampling FIOs In other words, we know an a priori bound of the smallest detail 𝑔 has. How to sample 𝐵𝑔? We need to fit WF 𝐵𝑔 in the smallest product Ω′ 𝐿 with 𝐿 a box or a parallelogram or a domain satisfying the tiling property w.r.t. some lattice. Then we sample on the reciprocal one. Now, if 𝐵 is a semi‐classical FIO or if WF 𝐵𝑔 was the classical WF set, then WF 𝐵𝑔 is related by WF 𝑔 by the canonical relation 𝐷 of 𝐵. It turns out that, aside from the zero section, this is also true in this case WF 𝐵𝑔 ∖ 0 ⊂ 𝐷 ∘ WF 𝑔 ∖ 0. Then the sampling of 𝐵𝑔 is determined by the geometry of C ∘ Ω 𝜊 𝐶 . (i) How to sample 𝐵𝑔 when 𝐵 is an FIO and WF 𝑔 ⊂ Ω 𝜊 𝐶 (for example)?

(ii) Resolution limit on 𝒈, given the sampling rate of 𝑩𝒈. (ii) Let us say we sample 𝐵𝑔 at same (limited) rate. This is what is done in

• practice. What is the smallest detail of 𝑔 we can possibly recover?

The way this is asked, it depends on whether 𝐵 is invertible, microlocally or not, in the first place. Let us say we sample 𝐵𝑔 on Ω′ at a rate 𝑡ℎ. Then if WF 𝑔 ⊂ C ∘ Ω 𝐶′, 𝐶′ we captured all “detail” of 𝐵𝑔 and did not lose any resolution about 𝑔 is 𝐵 is microlocally invertible in the first place. If 𝐵 is elliptic, then we can actually recover all the detail of 𝑔. If the above inclusion fails, there is aliasing. Note that the set on the right is not necessarily a product. This means resolution changing with location and direction.

(iii) Aliasing (iii) Let us say that 𝐵𝑔 is undersampled. How would that affect 𝑔? From now on, we assume a reconstruction with 𝜓̂𝜊 smooth (so it is a s.c. symbol); then 𝜓 ∈ 𝒯. An inspection of the reconstruction formula says that instead of 𝑔

, we get

𝐻𝑔

𝐻𝑔

• ∈𝒂

, with 𝐻 a semi‐classical FIO with a canonical relation 𝑇 ∶ 𝑦, 𝜊 ↦ 𝑦, 𝜊 2𝜌𝑙/𝑡 (and amplitude 𝜓̂𝑡𝜊/𝜌 when the phase is 𝜚 𝑦 𝑧 ⋅ 𝜊 2𝜌𝑗𝑙 ⋅ 𝑧/𝑡). The aliasing leaves 𝑦 in place but shifts the frequencies around (and changes their directions) as we saw above.

(iii) Aliasing This characterizes aliasing as an h‐FIO. The question however is what does aliasing of 𝐵𝑔 do to 𝑔 (not 𝐵𝑔)? Let us say that 𝐵 is elliptic and associated to a local diffeo. Then aliasing of the data gives us the following reconstruction 𝐵𝐻𝐵𝑔

• instead of 𝑔 (if there is no aliasing, only 𝑙 0 is there and 𝐻 Id). The

result is a sum of h‐FIOs with canonical relations 𝐷 ∘ 𝑇 ∘ 𝐷. 𝑇 are frequency shifts but the composition is not, in general! We can expect shifts in the 𝑦 variable as well!

(iv) Averaged measurements (iv) What if we sample 𝜚 ∗ 𝐵𝑔

with 𝜚𝑦 𝜗𝜚𝑦/𝜗?

Let’s say that 𝐵 is associated with a local diffeo 𝐷. The operator 𝑕 ↦ 𝜚 ∗ 𝑕 is an h‐ΨDO with a principal symbol 𝜚 . Then by Egorov’s theorem, 𝜚 ∗ 𝐵𝑔

𝐵𝑄𝑔 𝑃 ℎ 𝑔 ,

where 𝑄 is an h‐ΨDO with a principal symbol 𝜚 ∘ 𝐷. Depending on 𝜚 , this could blur 𝑔, so we are back to the problem (i) but we recover a blurred version 𝑄𝑔

• f 𝑔

.

This works if 𝜚 changes from a sampling point to another; then the samples became samples of 𝑅𝐵𝑔

with 𝑅 an h‐PDO. We can still apply Egorov’s

theorem. Finally, we can solve the inverse problem within the inverse problem: choose 𝑅 so that 𝑄 becomes what we want, say convolution with a fixed kernel.

Example: The Radon transform in “parallel geometry” We present some applications. Let ℛ be the weighted Radon transform in 2D ℛ𝑔 𝜕, 𝑞 𝜆 𝑦, 𝜕 𝑔 𝑦 𝑒𝑇.

⋅

Write 𝜕 𝜒 cos 𝜒, sin 𝜒. It is well known that ℛ is an FIO with canonical relation 𝐷 𝐷 ∪ 𝐷, where 𝐷 𝑦, 𝜊 arg 𝜊

• , 𝑦 ⋅ 𝜊/|𝜊|
• , 𝑦 ⋅ 𝜊
• , 𝜊
• .

Each 𝐷 is a diffeo, with inverse 𝑦, 𝜊 𝐷

𝜒, 𝑞, 𝜒

, 𝑞̂ given by 𝑦 𝑞𝜕 𝜒 𝜒 /𝑞̂𝜕 𝜒 , 𝜊 𝑞̂ 𝜕 𝜒 .

Radon transform: (i) Sampling Sampling ℛ (with 𝜆 1) was studied by BRACEWELL, RATTEY & LINDGREN and NATTERER. The proofs are not complete because they require two‐parameter asymptotics of Bessel functions. The results are actually asymptotic even though they are not stated as such. We will prove them for the weighted ℛ. Assume WF 𝑔 ⊂ 𝑦 𝑆, 𝜊 𝐶. Take 𝐷 of that and project it to the 𝜒 , 𝑞̂

• variables. We get

This is an elementary calculation knowing 𝐷. It is the same set obtained by the authors above. Σℛ𝑔

The smallest bounding box is 𝑆𝐶, 𝑆𝐶 𝐶, 𝐶. This determines relative sampling rates (before multiplying by ℎ) 𝑡 𝜌 𝑆𝐶 , 𝑡 𝜌 𝐶 To get a more efficient sampling set, notice that the plane can be tiled by translations by 𝑆𝐶, 𝐶 and 0,2𝐶 . Thus we can sample on the parallelogram grid 𝑡ℎ𝑋𝒂 with any 𝑡 1 and 𝑋 𝜌 𝑆𝐶 2 1 𝑆 . Those are 4/𝜌 times the phase volume x2. Ideally, it would be equal to it. The sampling points are 8/𝜌 times twice the phase volume of 𝑔 (since 𝐷 is a double cover, “twice” is the right target). Radon transform: (i) Sampling

Remember, Σ ℛ𝑔 is just a projection of WF ℛ𝑔 . We can do better than

• this. Here is ℛ𝑔 (with 𝜆 0) for 𝑔 consisting of 6 small Gaussians. Then we

superimpose WFℛ𝑔 (the projection) on it over each point. 𝜒 𝑞 We see that the cone shape of Σℛ𝑔 above is the union of all those cones, and it is also what happens at 𝑞 0. If we sum up (integrate) the areas of those “cones”, we get the sharp sampling count = 2x that for 𝑔. To get close to that, we need to split the rectangle into small horizontal strips and sample differently in each one of them. Radon transform: (i) Sampling

Radon transform: (ii) Resolution Let 𝑡 and 𝑡 be fixed. How would that affect 𝑔? From 𝐷, we have 𝑦 ⋅ 𝜊 𝜌/𝑡, 𝜊 𝜌/𝑡. The first one tells us that 𝑡_𝜒 does not restrict resolution much in the center (𝑦 0) or when 𝜊 ∥ 𝑦. On the other hand, 𝑡 restricts resolution uniformly (which is consistent with the intertwining property of 𝑒

and Δ).

• riginal

Undersampled in 𝑞 and reconstructed Undersampled in 𝜒 and reconstructed uniform blur plus non‐local artifacts best resolution Circular lines well

• resolved. Radial not.

Radon transform: (iii) Aliasing There were a lot of artifacts in those images due to aliasing of the data. As we discovered, the reconstructed image is an h‐FIO of the original (a sum of such). Assume we undersample angularly. The canonical relation on the data in this case is 𝑇 ∶ 𝑞, 𝜒 ↦ 𝑞̂, 𝜒 2𝜌𝑙/𝑡. No aliasing corresponds to 𝑙 0; and typically, in addition, we will see 𝑙 1 only. This affects 𝑔 as a sum of h‐FIOs with canonical relations 𝑦, 𝜊 ↦ 𝐷

∘ 𝑇 ∘ 𝐷 𝑦, 𝜊

𝑦 ≔ 𝑦 2𝜌𝑙 𝑡 𝜊 𝜊 , 𝜊 when the projection of 𝑦 to 𝜊 is in 𝜌/𝑡 , 𝜌/𝑡. So we see that 𝑦 shifts along 𝜊 but 𝜊 does not! This is in some sense the opposite to what happens in classical aliasing.

• riginal

reconstruction ℛ𝑔 undersampled (aliased) ℛ𝑔 In other words, the data ℛ𝑔 is aliased as usual: the effect is local and the frequencies are shifted. The reconstructed 𝑔 is affected differently: the frequencies are preserved but the phantom is shifted! and of the aliased ℛ𝑔. The Fourier transform of ℛ𝑔 The frequencies of ℛ𝑔 shift left and right. Radon transform: (iii) Aliasing 𝑞 𝜒

Radon transform: (iii) Aliasing Assume we undersample in the 𝑞 variable . A similar computations shows that both 𝑦 and 𝜊 shift.

• riginal

reconstruction ℛ𝑔 ℛ𝑔 undersampled in 𝑞 The higher frequency pattern disappears! It actually shifts out of the computational domain.

Radon transform: (iii) Averaged data Blur the data: 𝑅ℛ𝑔 with 𝑅 an h‐ΨDO with a compact support and principal symbol 𝑟 𝑦, 𝜊 . Egorov’s thm: 𝑔 is blurred by 𝑞 𝑦, 𝜊

• 𝑟 ∘ 𝐷
• 𝑟 ∘ 𝐷.

If 𝑅 is a convolution, i.e., 𝑟 𝜔𝑏 𝜒 | 𝑐 𝑞̂|, then 𝑞 𝑦, 𝜊 𝜔 𝑏 𝜊 𝑐 𝑦 ⋅ 𝜊 . The blur is uniform iff 𝑐 0, i.e., we average w.r.t. 𝑞 only. Well known. If 𝑏 0, the blur is position and direction dependent.

• riginal

ℛ𝑔 averaged in 𝑞 ℛ𝑔 averaged in 𝜒

Example: The Radon transform in “fan‐beam geometry” Lines parameterized by initial point 𝑆𝜕𝛽

• n the circle 𝑦 𝑆 and direction making

angle 𝛾 with the radial ray. This is diffeomorphic to the parallel geometry case. One can compute 𝐷 𝐷 ∪ 𝐷 as before. Assume as before WF 𝑔 ⊂ 𝑦 𝑆, 𝜊 𝐶. Then Σ ℛ𝑔 for all such 𝑔 can be computed as before. There is a new twist here: 𝐷 𝑦, 𝜊 may be aliased, 𝐷 𝑦, 𝜊 may be not! 𝑔 ℛ𝑔 |ℛ𝑔 |

Fan‐beam geometry: (i) Sampling Σℛ𝑔 Example: 𝑔 Computed |ℛ𝑔 | The double rectangle above determines the sampling requirements on a rectangular grid: 𝑡 𝜌 𝑆𝐶 , 𝑡 𝜌 2𝑆𝐶 . This is 𝜌 the parallel case (on a rectangular grid), and 8 more than the Weyl optimum. We can tile the plane with translations by 𝑆𝐶, 0 and 0,2𝑆𝐶 as we showed

• above. This doubles the step sizes. We still have twice the optimum number.

Fan‐beam geometry: (i) Sampling A better representation of WFℛ𝑔: The figure we got earlier is the “worst scenario case”: the union of all those. The greatest sampling requirements are near 𝛾 0. One can sample over small parallel strips with different steps. This would bring us asymptotically close to the Weyl minimum. ℛ𝑔 with WF ℛ𝑔 |, at a few points. Σℛ𝑔

𝛽 𝛾

Fan‐beam geometry: (ii) Resolution Assume we sample ℛ𝑔 with relative rates 𝑡 and 𝑡. What does this imply for 𝑔? We need to understand 𝐷 𝛽, 𝛾 ∈ 0,2𝜌

• ,
• ; 𝛽
• , 𝛾
• .

The red bars represent the frequency cutoff caused by 𝑡. The black circles correspond to 𝑡. The frequencies in the intersection of the circles are the non‐aliased ones. The reason: 𝐷 is 2‐to‐1. As a result: there is more resolution in the data, if extracted properly! The diagram is rotationally symmetric.

Fan‐beam geometry: (iii) Aliasing If the sampling rates 𝑡 and 𝑡 are not small enough, aliasing occurs a

• before. A direct computation shows that 𝑦 shifts along the direction of 𝜊

and 𝜊 changes magnitude but not direction.

Fan‐beam geometry: (iv) averaged data Blur the data: 𝑅ℛ𝑔 with 𝑅 an h‐ΨDO with a compact support and principal symbol 𝑟 𝑦, 𝜊 . Egorov’s thm: 𝑔 is blurred by 𝑞 𝑦, 𝜊

• 𝑟 ∘ 𝐷
• 𝑟 ∘ 𝐷.

If 𝑅 is a convolution, i.e., 𝑟 𝜔𝑏 𝜒 | 𝑐 𝑞̂|, then 𝑞 𝑦, 𝜊 1 2 𝜔 𝑏 𝑦 ⋅ 𝜊 𝑐 𝑦 ⋅ 𝜊 𝑆 𝜊 𝑦 ⋅ 𝜊 1 2 𝜔 𝑏 𝑦 ⋅ 𝜊 𝑐 𝑦 ⋅ 𝜊 𝑆 𝜊 𝑦 ⋅ 𝜊 . The levels sets are as in the resolution analysis. Unlike in the parallel geometry case, those terms are different! They blur differently and then the results are added.

Fan‐beam geometry: (iv) averaged data (a) original (b) 𝑔 with ℛ𝑔 averaged in α (c) 𝑔 with ℛ𝑔 averaged in 𝛾 (b) Radial blur. Corresponds to the red bars. (c) Sum of two blurs. Corresponds to the circles. (d) Source restricted to the right semi‐circle. We get better resolution near the source (direction dependent). Corresponds to one of the circles

• nly.

The resolution diagram (d) 𝑔 with partial ℛ𝑔 averaged in 𝛾

Fan‐beam geometry: (iv) averaged data To understand better the effect we see in (c) above, here is another test function: a 5 5 array of almost delta‐like Gaussians. Reconstructed 𝑔 with ℛ𝑔 averaged in 𝛾, i.e., data: 𝜔 𝛾 ∗ ℛ𝑔 A magnified crop A density plot of the computed symbol at 0.88𝑆, 0 𝑞 𝑦, 𝜊 1 2 𝜔 𝑦 ⋅ 𝜊 𝑆 𝜊 𝑦 ⋅ 𝜊 1 2 𝜔 𝑦 ⋅ 𝜊 𝑆 𝜊 𝑦 ⋅ 𝜊

Fan‐beam geometry: anti‐aliasing Bluring ℛ𝑔 can be used to avoid aliasing.

• riginal

𝑔 with ℛ𝑔 undersampled in 𝛾 𝑔 with ℛ𝑔 blurred in 𝛾, then sampled A lot of aliasing artifacts spreading everywhere The anti‐aliasing removes most of the artifacts but softens the image a bit more.

Thermo‐acoustic Tomography Model: 𝑣 𝑑 𝑦 Δu 0 in 𝐒 𝐒, 𝑣| 𝑔 𝑦 , 𝑣| 0, where 𝑣𝑢, 𝑦 is the sound pressure, 𝑑 𝑦 is the sound speed equal to 1

• utside Ω.

Ω 𝑢 𝑈 We measure Λ𝑔 𝑣| , . There is no actual boundary, the waves propagate to the whole space. Clearly, Λ is a classical FIO. Another studied model is when there are Neumann b.c. and the waves reflect. Λ is still a classical FIO. 𝑣𝑢, 𝑦

Thermo‐acoustic Tomography Canonical relation 𝐷 𝐷 ∪ 𝐷. 𝐷𝑦, 𝜊 𝐷𝑦, 𝜊 For every 𝑢, 𝑦 ∈ 0, 𝑈 𝜖Ω, the range of the possible 𝐷 𝑦, 𝜊 (when 𝜊 𝐶) is the characteristic cone projected there. Their union gives ΣΛ𝑔, which is the cone above. Here, 𝑁 max 𝑑 because we have 𝑁 1 if we use the metric 𝑑𝑒𝑦 on 𝑈∗Ω but we have 𝜊 𝐶 with 𝜊 Euclidean (as determined by the sampling rate). The induced metric on 𝜖Ω would results in another constant 𝑁 but when 𝜖Ω is piecewise flat, 𝑁 1. 𝑦 𝜊 ΣΛ𝑔

Thermo‐acoustic Tomography It turns out that the geometry plays no role! Only max 𝑑 and the shape

• f 𝜖Ω matter!

𝑔 Fast region inside, no caustics (before reflections). Much higher frequencies, needs finer sampling. Slow region inside, caustics. Waves focus

• n 𝜖Ω. Low sampling/ discretization

requirements. Compute Λ𝑔in the model with reflections. The r.h.s. shown. 𝑢 𝑦 𝑢 𝑦