Back to Teaching page

UW RTG IPDE Summer School 2013 (website)

Labs 1(07/09), 2(07/11) and 3(07/16): See Peter Caday's webpage.

Lab 4 (Th 07/18): Ray transforms: different families of curves, fun with FIO's.

CODES
In this session, we will not construct an exact inversion. We will simply backproject the data for different types of geometries of curves and see what happpens. We will see how the choice of integration curves is the only thing that matters in the transport of singularities (and not the integration weights as long as they are smooth and non-zero), and how to typically compute these singularities. Moreover, we will see how errors on the geometric features of the operator (in this case, the family of lines, in the TAT case, the sound speed) generates artifacts that, in some cases, can be analytically computed.
  1. Construct a family of curves whose geometry depends on one parameter $R$. We will construct in class a family that varies continuously from straight lines to the lines of the Poincare unit disk of constant negative curvature.
  2. Fill out the UWgeodesic and UWhyperS functions. The first one returns a few points describing the unique curve in the $R$-hyperbolic model of "coordinates" $(s,\theta)$, i.e. the unique curve $\gamma_{s,\theta}(t)$ satisfying $\gamma_{s,\theta}(0) = s\vec{\theta}$ and $\dot\gamma_{s,\theta} (0) = \vec{\theta}^\perp$. The second function is necessary for backprojection: for a given point $x$ and a direction $\theta$, it returns the unique $s$ such that the geodesic of "coordinates" $(s,\theta)$ passes through $x$. Knowledge of the function $s(x,\theta)$ is enough to implement the backprojection operators. Run Hyper1 as main program to test the exactness of your formulas.
  3. Use these functions to run the first part of the Hyper2 script and notice the difference between regular and "hyperbolic" Radon Transform (pick a value of $R$ between $1$ and $1.5$ to notice the difference). You will also need the following scripts: hyperRadon, hyperbackproj, laplacian (one applies the Laplace operator to the backprojected data once to sharpen the singularities).
  4. Intertwining operators with different geometries. Play with the codes (second part of the Hyper2 script), see how singularities are transformed and notice the artistic potential of microlocal analysis.
  5. Open questions:
    1. Design another family of curves such that $s(x,\theta)$ is analytically defined, code geodesic and hyperS functions accordingly and repeat the above questions (suggestions: pertubation of straight lines, parabolas).
    2. For the example above, compute the graph relation of the compositions of FIOs, obtained after solving \[ s_1(x,\theta) = s_2 (y,\theta)\quad \text{and} \quad \partial_\theta s_1 (x,\theta) = \partial_\theta s_2 (y,\theta). \]
    3. "Calibration problem": Suppose one computes the ray transform with one given family of curves and the backprojection with another. Could one infer something about how far off the families of curves are from knowledge of pairs of initial images and their corresponding reconstructed images ?
Additional functions:

Lab 5 (Tue 07/23): Thermoacoustic Tomography: More building blocks toward the inversion formula.

CODES
In this session, we will cover the main two steps of the approximate inversion procedure: (i) constructing a harmonic continuation of the boundary data at time $T$ and (ii) solving a backward wave equation with boundary conditions. We will also see how to solve eikonal equations, which in turn allows us to compute the minimum time for injectivity of the TAT problem. Solutions:

Lab 6 (Th 07/25): Thermoacoustic Tomography: Implementation of the Neumann series.

CODES
Now that all the building blocks are constructed, we are now ready to implement the Neuman series, which iteratively improves the reconstruction. Recall that the reconstruction formula reads \[ f - Kf = A\Lambda f,\] where $\Lambda f$ denotes the data or forward map, $A$ denotes approximate inversion (comprising of a harmonic extension and the resolution of backward wave problem) and $K$ denotes the "error" operator (which in fact can be defined as $K := Id - A\Lambda$). When the observation time $T$ is large enough, the operator $K$ is proved theoretically to be a contraction, so that the relation above may be inverted via the following Newmann series \[ f = (Id-K)^{-1} A\Lambda f = \sum_{k=0}^\infty K^k A\Lambda f = \sum_{k=0}^\infty (I-A\Lambda)^k A\Lambda f, \] and the purpose of today's session is to implement a finite number of terms of this series.