Machine-Learning-Homework 7 Solved

Setting
You nd yourself in a two-dimensional world in the role of a medical doctor. A patient (we will call him H.S. here for anonymity) walks into your hospital and complains about a headache. A close examination of the head provides no exterior clues as to the cause of the pain, and you soon realize that you will have to perform computer tomography. Unfortunately, the computer for processing recorded data is broken. Nevertheless, you let H.S. enter the tomograph and take X-ray images of his head from several angles. You send H.S. back to the waiting room and start to think about how you might process the scans using only your private laptop.

Preliminaries
For the scans, H.S.’s 2D head was placed between an X-ray source and a 1D sensor array. The sensor elements record the accumulated signals of many parallel rays passing through the head (see left side of gure 1). To reconstruct all the internal structures, multiple measurements need to be taken from di erent angles α. The resulting tomogram is an image µ ∈ RM×M where each pixel encodes the local X-ray absorption coe cient of the patient’s head (bright for bones, gray for soft tissue, black for air). The least-squares reconstruction algorithm will create a attened version β ∈ RD of µ that needs to be reshaped into a 2D image to be displayed. The pixels µja,jb are related to the vector elements βj by the formula j = ja + Mjb, as shown in gure 1.

The sensor consists of a 1-dimensional array of Np sensor elements that register parallel X-rays, and the emitter-detector setup is rotated into No di erent orientations αio. The measured intensities are commonly displayed as a 2D image called a sinogram, see gure 3. For least-squares reconstruction, we atten the sinogram s ∈ RNo×Np into the response vector y ∈ RN. Its elements yi are de ned by the relation i = ip + Npio, where ip and io are the indices of the sensor elements and sensor orientations respectively.

We can now describe the relation between the tomogram β and the measured intensities y as a linear projection with projection matrix X ∈ RN×D:

X · β = y (1) Intuitively, entry Xij should be interpreted as a weight which encodes how much the sensor response yi is in uenced by the absorption at pixel βj.

In the continuous domain, the measured intensity Isensor of a single ray r depends on how much of the original radiation Isource is absobed by the material between emitter and sensor. This can be modelled as

               Isensor = Isource  ,                                                                                          (2)

where µ describes the material’s absorption properties at each spatial position (a,b) and can be seen as a continuous version of the discrete tomography image µ de ned above. For simplicity, we assume that the raw sensor measurements are already subjected to preprocessing which computes the logarithm of this formula, so that the response is de ned by

 

 source

                     y = ln                                                                                                                                              (3)

sensor

Z =
µ(a,b)dadb.
(4)
a,b∈r

For practical computation, we need to discretize this formula. Each sensor element is essentially a bin collecting radiation from many parallel rays. We assume that a ray intersecting the sensor array inbetween two bins distributes its intensity among both. The ratio of distribution is determined in a linear fashion by how close the intersection is to the center of either bin.

Example: Consider the ray passing through β19 at angle αio, as illustrated in gure 1, which intersects the sensor array at position 4.35 (between elements 4 and 5). Sensor element 4 captures 65% of this ray’s intensity and sensor element 5 captures 35%. This is encoded in the weight matrix X by setting X4+Npio,19 = 0.65 and X5+Npio,19 = 0.35.

Thanks to this de nition, the integral in equation 4 can be discretized into a sum, and writing the sum in matrix notation gives the linear system of equation 1.

 

Figure 1: Sketch of the projection of a 5×5-image µ onto a 7 pixel sensor array that has been rotated around the center of the image by the angle αio (left). Small numbers re ect the indices j within β and ip within one slice of s, respectively. While both the X-ray source and the sensors are naturally situated outside the tissue, the projection can be made easier by virtually placing the sensor at the image center and considering the rays to come from both directions (right). In any case, for each individual coordinate in β, we trace the ray passing through that point and calculate how much of its intensity is collected by which sensor pixels.

1           Constructing the matrix X 
With only one projection, solving equation 1 is ill-posed and a useful reconstruction of β is impossible. The problem becomes well-posed by using su ciently many projections of β under di erent angles α, whose recordings are just concatenated in the attened measurement vector y or the 2D sinogram ( gure 3), respectively. Constructing the corresponding weight matrix X is the main di culty of this exercise.

Implement a function X = construct_X(M, alphas, Np = None) which, given the desired tomogram size D = M × M, a list of measurement angles (α0,α1,...) (in degrees) and an optional sensor resolution Np, returns the matrix X, which completely describes the setup of our simpli ed CT scanner.

If Np is not given, choose a value large enough to t the diagonal of the image β, like . Choosing an odd number of sensor elements will make it easier to align the image coordinates with the coordinates of the sensor array: You can de ne a common coordinate origin and then nd ray intersections by projecting the tomogram’s pixel coordinates along the current ray orientation onto the rotated sensor array, see gure 1.

In the examination scenario for patient H.S., we have M = 195, i.e. the tomogram has D = 38025 pixels. The sensor has Np = 275 bins, and intensities have been measured at No = 179 projection angles, resulting in a response vector of size N = 49225. Thus, matrix X has shape 49225 × 38025 (almost 2 billion entries), and you will run into trouble if you implement construct_X() naively.

Fortunately, the intensity of each ray is only distributed among the two sensor elements closest to the ray’s intersection with the sensor. Consequently, each image pixel can contribute to no more than two sensor elements under any given measurement angle (see right side of gure 1) and the vast majority of entries of X are actually zero. This is called a sparse matrix, and the scipy module scipy.sparse[1] provides powerful tools to take advantage of this property: it allows you to save memory by o ering spare matrix classes that store only the non-zero entries, and to skip unnecessary calculations by o ering specialized implementations of linear algebra functionality. Speci cally, the class coo_matrix ( COOrdinate format) is recommended to create X, and the class csc_matrix ( Compressed Sparse Column format) is best for solving the linear system. Converting between di erent matrix types is easy and results in faster code than using the same type throughout.

However, using a sparse matrix is not enough to make your code e cient. You should also vectorize the function construct_X() (cf. exercise 1). A good way of doing this is to create an array C ∈ R2×D holding the coordinates of the tomogram’s pixel centers, i.e.

C0j = a0 + ja · h (the a coordinate of pixel j)

C1j = b0 + jb · h (the b coordinate of pixel j)

with h being the pixel distance and j = ja + Mjb as de ned earlier (numpy.mgrid is useful here). In this case, the pixel distance for both image and sensor is simply h = 1. The current orientation can be expressed by a unit vector n along the rotated sensor, and the projection p of each pixel onto the sensor is simply a linear projection according to

p = n>C + s0,

where s0 is the distance between the rst sensor element and the sensor’s coordinate origin. The weights X for each angle can now be computed from p by numpy’s vector functions. In the end, the valid elements of X are represented by their indices (i,j) and the corresponding weights. The constructor of class coo_matrix expects this information to be provided by three 1-dimensional arrays i_indices, j_indices, and weights, resulting in the call

X = coo_matrix((weights, (i_indices, j_indices)), shape=(N, D), dtype = numpy.float32),

 

Figure 2: Visualization of matrix X for M = 10, Np = 15 and three projections at angles (−33,1,42). You can also nd this matrix as ‘hs_tomography/X_example.npy’ on MaMpf.

 

Figure 3: Sinogram visualization of y. Each column shows the sensor response

for a di erent angle .

where we used element type float32 to save even more memory. Check the scipy documentation for more details about how this works. In our example implementation, constructing the sparse 49225 × 38025 matrix takes about 5 seconds on a standard laptop.

As in real-world scenarios, there is some ambiguity left to this task: Where is the origin of the coordinate system, which direction is α = 0, and how is y oriented? The answers to these questions determine the appropriate indexing order and sign of the expressions, and guring out these details is part of your task (emulating what often happens in the real world). To check whether your matrix construction is correct, visualize the result of construct_X(10, [-33, 1, 42]) with pyplot.imshow() and compare it to the matrix shown in gure 2, which you can download from MaMpf as ‘hs_tomography/X_example.npy’[2]. Note that you will have to convert your sparse matrix to a dense numpy-array for visualization.

2           Recovering the image 
The material for this sheet contains the list of angles and the measured sensor data y for two versions of the experiment. The smaller one was created with M = 77, Np = 109 and 90 projection angles. The larger one was created with M = 195, Np = 275 and 179 projection angles. The exact set of angles and the corresponding response vectors y can be downloaded from MaMpf in the les ‘alphas_77.npy’ and ‘y_195.npy’ resp. ‘alphas_195.npy’ and ‘y_77.npy’ (in folder ‘hs_tomography’). We recommend the smaller version of the data for debugging your code. However, its resolution is insu cient to diagnose the cause for your patient’s headache. If you can recognize the head in the small tomogram, your code is probably correct, and you should switch to the larger version.

Exactly how many non-zero entries does X have? Use scipy’s tools to nd out and report the sparsity of X.

With the ability to construct the matrix X and having obtained a set of projections y, you are now able to reconstruct the original image β. In theory, you could try to solve the equation system directly via the pseudo-inverse

β = (X>X)−1 · X> · y,

but this is not recommended because it ignores sparsity and is therefore very slow.

Fortunately, scipy already gives you access to e cient solvers for sparse linear equation systems! Find out how to use scipy.sparse.linalg.lsqr() to obtain the least-squares solution to your problem. The solver’s low default tolerance parameters atol = 1e-08 and btol = 1e-08 lead to high quality solutions, but also long computation times on the large version of the data. A good trade-o might be atol = btol = 1e-05, but you can reduce the tolerance once your are con dent that your code is correct.

Reconstruct the tomogram and plot it as a 2D image. Give a diagnosis on what causes H.S.’s headache and propose a treatment.

If you didn’t manage to correctly construct X in task 1, you can nd a precomputed weight matrix for the smaller version of the experiment in ‘hs_tomography/X_77.npy’ and use it for this task. You can load it with X = np.load(’hs_tomography/X_77.npy’, allow_pickle=True) and convert it to a sparse matrix with X = scipy.sparse.csc_matrix(X.all()).

3           Minimizing the radiation dose
As a doctor, you do not want to expose your patients to unnecessary radiation. Try to reduce the number of projection angles in a sensible way and visualize how this changes the quality of the reconstruction. In the case of H.S., what would you say is the minimal number of projections that still allows you to resolve the cause of his headache?


 
[1] Install scipy via conda as usual.
[2] 2‘.npy’ is numpy’s matrix                                   le format. You can load it with X_ex = np.load(’hs_tomography/X_example.npy’)