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.
- 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.
- 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.
- 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).
- 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.
- Open questions:
- 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).
- 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). \]
- "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.
- Harmonic extension:
- Complete the diffOps function, whose aim is to construct a discretized version of the Laplacian operator with Dirichlet boundary conditions, after constructing a discretized gradient $(\partial_x,\partial_y)$.
- Test the correctness of your code by running demoPoissonSolver and checking that the solution of your Poisson equation matches that of some problem whose solution is known exactly.
- The backpropagation problem:
- Complete the backwave function, whose aim is to solve a backward wave problem on a bounded domain with prescribed boundary conditions and final conditions. The final condition $u|_{t=T}$ is the input uend while the initial velocity satisfies $\partial_t u|_{t=T}=0$.
- Check the exactness of your code by running the DemoSingleInversion script, which computes a full (forward solver)-(harmonic extension)-(backward solver) sequence. With the parameter values already plugged in for the sound speed, initial condition, observation time, etc... the reconstructed function should be visually close to the initial one.
- Experiments: When the functions above are coded and exact, you can now run DemoSingleInversion with different values of the parameters and verify some principles that were predicted theoretically.
- What happens when $T$ is too small for injectivity ?
- What happens when $T$ is large enough for injectivity but too small for stability ?
- What happens when the sound speeds are chosen to be different in the forward solver and the backward solver ?
- Compare the quality of reconstructions for a smooth input versus a discontinuous input, everything else being equal.
- Compare this method with the time-reversal method by setting postproc = 'TR' in the script (this causes the script to NOT carry out the harmonic extension step).
- Compare the impact of trapping/non-trapping speeds on reconstructions.
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.
- Using the functions wavePML and backwave, complete the script DemoIterative in order to implement the Neumann series. At this point, testing that the Neumann series does a good job consists in checking that the iterations improve if the observation time has been chosen large enough. One may check convergence by visualizing the $L^2$ norm of the error as iterations increase.
- For a given sound speed map, use the eikonal function to get an estimate of the injectivity time.
- Carry out some of the experiments listed in the Lab 5.