ScanPyramid Midterm Project.

Background.

In this project, we are modeling the project of muon imaging of the great pyramid. In muon imaging, detectors are placed around the structure to be investigated. Muons arrive at random angles as cosmic rays. However, when the muons pass through dense substances (such as rock), some of them are absorbed or deflected. 

Experimenters can record the number of muons that arrive on a given flight line and use these counts to gradually build up an integral of the density of the structure along the flight line. Given enough of these flight lines, the internal structure of the object can be recomputed.

The data.

Since this is an exercise, we’re going to use a two-dimensional model of a pyramid instead of a 3d model. The density of the pyramid is a function μ(x,y) where x ∈ [-1500,1499] and y ∈ [-1273,1273]. This rectangle is divided into an array of 3000 x 2547 square pixels, where μ(x,y) = 1 if that square is occupied by solid material and μ(x,y) = 0 if the corresponding square is empty.

We will parametrize the space of possible flight lines in the plane by two coordinates θ and d, where the line (θ,d) has equation d = x cos θ + y sin θ. The lines will have θ values between -π/2 and π/2, and d values between plus or minus

Sqrt[3000^2 + 2547^2]/2

N[((3 (Sqrt[1720801]))/2)]

which is the diagonal of the image (in pixels). The coordinate system is centered at the center of the image. In these coordinates, the pyramid defines a density function μ(x,y) which is 1 for pixels within the solid portion of the pyramid and 0 for empty pixels.

For a selection of lines with uniformly random θ and d coordinates, we’ve computed the integral

R(θ,d) = (∫_(-∞)^∞)μ(-t sin θ + d cosθ,t cos θ + d sin θ)dt

that is, the total density along the line (θ,d).

The data has been normalized so that the largest value of the integral is 1 and the smallest is 0. The file containing the data is located on the course webpage, and called Pyramid-Scan-Data.csv. 

It contains a gracious plenty of data (half a million samples), so you probably want to test your methods on a smaller sample of the data before running on the entire dataset. As is so often the case with data, the columns in the file are unlabeled. You’ll have to look at the values in each column to decided which is which.

Mission.

Your mission is to recover the density function μ(x,y) as well as you can from the information in Pyramid-Scan-Data using one of the methods listed below.

Full Credit Methods.

There are a number of known methods for recovering the density map μ(x,y) from the R(θ,d) function. You must implement one of these using numerical integration and differentiation code of your own devising, based on the methods we’ve learned in class.

The Direct Formula.

μ(x,y) = ((-1)/(2(π^2)))(∫_(-π/2)^(π/2))(∫_(-∞)^∞)((∂R(θ,d))/(∂d))(1/(d - x cos θ - y sin θ))dddθ

Notice that the function R(θ,d) is 0 (and hence has partial derivative 0) for d values larger and smaller than half the diagonal of the original image, so the limits of integration on the inner integral amount to “integrate over all the nonzero values”.

The Fourier Method.

Compute the Fourier transform 

(OverHat[R])(θ,d) = (∫_(-∞)^∞)R(θ,x) (e^(2π i x d))dx

Suggested first step. First, isolate all the data corresponding to a single θ value (for instance, θ = 0, which refers to vertical lines). You should have enough different x values present to make a plot of the function R(0,x). What should this function look like? You should have enough data to construct a pretty good interpolation of R(0,x), and compute the Fourier transform of that function. Depending on how you interpolate, you may be able to get an explicit formula for the Fourier transform function (OverHat[R]). But you should definitely be able to get values at enough points to make a plot of (OverHat[R])(0,d).

Suggested second step. If you’re confident in your computation for (OverHat[R])(0,d), automate the work above and construct a table of data for (OverHat[R])(θ,d) for a table of θ values. This will give you a new dataset for the final integration. (It’s up to you and your computer exactly how big this dataset is. But I’d eventually make it pretty large.)

Define a new function on the u-v plane in polar coordinates, where u = d cos θ and v = d sin θ by

F(u,v) = F(d cosθ, d sin θ) = (Overscript[R, ⋀])(θ,d)

Then it turns out to be the case that

μ(x,y) = ∫∫F(u,v)(e^(2 π i (x u + y v)))dudv

where the integral is taken over everywhere that F(u,v) is nonzero.

The cutoff function method.

Amazingly, it’s also true that 

μ(x,y) = (lim_(z -> 0))(∫_(-π/2)^(π/2))(∫_(-∞)^∞)R(θ,t-x cos θ - y sin θ) Γ(z,t)𝕕 tdθ

where 

Γ(z,t) = (1/(π (z^2))), if -z ≤ t ≤ z

and 

Γ(z,t) = (1/(π (z^2)))(1 - (1/(Sqrt[1 - (z^2)/(t^2)]))) , for |t| > z.

Partial Credit Methods.

For partial credit, you can use Mathematica’s built-in integrators instead of writing your own. In partial, you can use Mathematica’s built-in Fourier transform methods to compute the Fourier transforms in method 2.

Implementation notes.

One of the largest issues in this kind of integration is the question of how to deal with missing data values. Of course, 1-d polynomial interpolation is your friend here, but you’re encouraged to read the Mathematica tutorial “Approximate Functions and Interpolation” and consider your other options as well. 

You’re going to have to synthesize data on a rectangular grid in θ, d coordinates in order to use any of our methods going forward to approximate the numerical derivatives and integrals. (Why?)

Start by recovering ONE pixel value of μ(x,y) that you can guess before leaping into computing millions of them.

The Piecewise function in Mathematica can be very helpful when you’re working with functions like Γ(z,t) which are defined, well, piecewise.

If your method is slow, it may be useful to know that the machines in the classrooms are pretty good and have Mathematica installed. Don’t burn them out; if they’re running, you need to be there.

Rules.

The rules are those of a term paper. This is not an exam. Your group may consult any inanimate source of information (that is, books, Mathematica documentation, Google searches, but not crowdsourcing sites such as StackExchange). No people, except for me.

Data will be posted by Midnight on Tuesday October 24. Completed Mathematica notebooks are due by email by Midnight on Tuesday November 7. 

One Mathematica notebook per group. 

Your notebook should discuss the numerical methods you tried and developed, document your code, cover the material in any reading or research you did on the problems involved, and finally present a single solution as a bitmapped image which is 3000 x 2547 pixels. You may also describe any areas inside the pyramid which you believe to contain voids even if you can’t reconstruct an entire image.
You are allowed to use the Mathematica image processing functions to improve your results. 

You may ask me questions about your term paper project after class and during office hours.

Acknowledgements.

Thanks to Blake Pinto for pointing out an error in the original dataset, to Project DANGEROUS WRENCH for pointing out that the formula for the Fourier transform was wrong. DANGEROUS WRENCH continued to help contribute edits and waypoints to the Fourier method.