Tomography

Tomography is the process of reconstructing a density distribution from given integrals over sections of the distribution. In our example, we will work with tomography on black and white images. Suppose $x$ be the vector of $n$ pixel densities, with $x_j$ denoting how white pixel $j$ is. Let $y$ be the vector of $m$ line integrals over the image, with $y_i$ denoting the integral for line $i$. We can define a matrix $A$ to describe the geometry of the lines. Entry $A_{ij}$ describes how much of pixel $j$ is intersected by line $i$. Assuming our measurements of the line integrals are perfect, we have the relationship that

$$y = Ax$$

However, anytime we have measurements, there are usually small errors that occur. Therefore it makes sense to try to minimize

$$\|y - Ax\|_2^2.$$

This is simply an unconstrained least squares problem; something we can readily solve!

In [2]:
using Convex, Gadfly, SCS

line_mat_x = readdlm("tux_sparse_x.txt")
line_mat_y = readdlm("tux_sparse_y.txt")
line_mat_val = readdlm("tux_sparse_val.txt")
line_vals = readdlm("tux_sparse_lines.txt")

# Form the sparse matrix from the data
# Image is 50 x 50
img_size = 50
# The number of pixels in the image
num_pixels = img_size * img_size

line_mat = spzeros(3300, num_pixels)

num_vals = length(line_mat_val)

for i in 1:num_vals
  x = int(line_mat_x[i])
  y = int(line_mat_y[i])
  line_mat[x + 1, y + 1] = line_mat_val[i]
end

pixel_colors = Variable(num_pixels)
# line_mat * pixel_colors should be close to the line_integral_values
# to reflect that, we minimize a norm
objective = sumsquares(line_mat * pixel_colors - line_vals)
problem = minimize(objective)
solve!(problem, SCSSolver(verbose=0))

rows = zeros(img_size*img_size)
cols = zeros(img_size*img_size)
for i = 1:img_size
  for j = 1:img_size
    rows[(i-1)*img_size + j] = i
    cols[(i-1)*img_size + j] = img_size + 1 - j
  end
end

# Plot the image using the pixel values obtained!
p = plot(
  x=rows, y=cols, color=reshape(evaluate(pixel_colors), img_size, img_size), Geom.rectbin,
  Scale.ContinuousColorScale(Scale.lab_gradient(color("black"), color("white")))
)
Out[2]:
x -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 -50 -48 -46 -44 -42 -40 -38 -36 -34 -32 -30 -28 -26 -24 -22 -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70 72 74 76 78 80 82 84 86 88 90 92 94 96 98 100 -50 0 50 100 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 y 0 300 100 -100 200 Color -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 -50 -48 -46 -44 -42 -40 -38 -36 -34 -32 -30 -28 -26 -24 -22 -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70 72 74 76 78 80 82 84 86 88 90 92 94 96 98 100 -50 0 50 100 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100