You can read an overview of this Numerical Linear Algebra course in this blog post. The course was originally taught in the University of San Francisco MS in Analytics graduate program. Course lecture videos are available on YouTube (note that the notebook numbers and video numbers do not line up, since some notebooks took longer than 1 video to cover).

4. Compressed Sensing of CT Scans with Robust Regression¶

The term broadcasting describes how arrays with different shapes are treated during arithmetic operations. The term broadcasting was first used by Numpy, although is now used in other libraries such as Tensorflow and Matlab; the rules can vary by library.

From the Numpy Documentation:

Broadcasting provides a means of vectorizing array operations so that looping
occurs in C instead of Python. It does this without making needless copies of data
and usually leads to efficient algorithm implementations.

The simplest example of broadcasting occurs when multiplying an array by a scalar.

In [718]:
a = np.array([1.0, 2.0, 3.0])
b = 2.0
a * b

Out[718]:
array([ 2.,  4.,  6.])
In [128]:
v=np.array([1,2,3])
print(v, v.shape)

[1 2 3] (3,)

In [129]:
m=np.array([v,v*2,v*3]); m, m.shape

Out[129]:
(array([[1, 2, 3],
[2, 4, 6],
[3, 6, 9]]), (3, 3))
In [133]:
n = np.array([m*1, m*5])

In [134]:
n

Out[134]:
array([[[ 1,  2,  3],
[ 2,  4,  6],
[ 3,  6,  9]],

[[ 5, 10, 15],
[10, 20, 30],
[15, 30, 45]]])
In [136]:
n.shape, m.shape

Out[136]:
((2, 3, 3), (3, 3))

In [48]:
m+v

Out[48]:
array([[ 2,  4,  6],
[ 3,  6,  9],
[ 4,  8, 12]])

Notice what happens if we transpose the array:

In [49]:
v1=np.expand_dims(v,-1); v1, v1.shape

Out[49]:
(array([[1],
[2],
[3]]), (3, 1))
In [50]:
m+v1

Out[50]:
array([[ 2,  3,  4],
[ 4,  6,  8],
[ 6,  9, 12]])

When operating on two arrays, NumPy compares their shapes element-wise. It starts with the trailing dimensions, and works its way forward. Two dimensions are compatible when

• they are equal, or
• one of them is 1

Arrays do not need to have the same number of dimensions. For example, if you have a $256 \times 256 \times 3$ array of RGB values, and you want to scale each color in the image by a different value, you can multiply the image by a one-dimensional array with 3 values. Lining up the sizes of the trailing axes of these arrays according to the broadcast rules, shows that they are compatible:

Image  (3d array): 256 x 256 x 3
Scale  (1d array):             3
Result (3d array): 256 x 256 x 3

Review¶

In [165]:
v = np.array([1,2,3,4])
m = np.array([v,v*2,v*3])
A = np.array([5*m, -1*m])

In [166]:
v.shape, m.shape, A.shape

Out[166]:
((4,), (3, 4), (2, 3, 4))

Will the following operations work?

In [159]:
A

Out[159]:
array([[[ 5, 10, 15],
[10, 20, 30],
[15, 30, 45]],

[[-1, -2, -3],
[-2, -4, -6],
[-3, -6, -9]]])
In [158]:
A + v

Out[158]:
array([[[ 6, 12, 18],
[11, 22, 33],
[16, 32, 48]],

[[ 0,  0,  0],
[-1, -2, -3],
[-2, -4, -6]]])
In [167]:
A

Out[167]:
array([[[  5,  10,  15,  20],
[ 10,  20,  30,  40],
[ 15,  30,  45,  60]],

[[ -1,  -2,  -3,  -4],
[ -2,  -4,  -6,  -8],
[ -3,  -6,  -9, -12]]])
In [168]:
A.T.shape

Out[168]:
(4, 3, 2)
In [169]:
A.T

Out[169]:
array([[[  5,  -1],
[ 10,  -2],
[ 15,  -3]],

[[ 10,  -2],
[ 20,  -4],
[ 30,  -6]],

[[ 15,  -3],
[ 30,  -6],
[ 45,  -9]],

[[ 20,  -4],
[ 40,  -8],
[ 60, -12]]])

Sparse Matrices (in Scipy)¶

A matrix with lots of zeros is called sparse (the opposite of sparse is dense). For sparse matrices, you can save a lot of memory by only storing the non-zero values.

Another example of a large, sparse matrix:

There are the most common sparse storage formats:

• coordinate-wise (scipy calls COO)
• compressed sparse row (CSR)
• compressed sparse column (CSC)

Let's walk through these examples

There are actually many more formats as well.

A class of matrices (e.g, diagonal) is generally called sparse if the number of non-zero elements is proportional to the number of rows (or columns) instead of being proportional to the product rows x columns.

Scipy Implementation

• To construct a matrix efficiently, use either dok_matrix or lil_matrix. The lil_matrix class supports basic slicing and fancy indexing with a similar syntax to NumPy arrays. As illustrated below, the COO format may also be used to efficiently construct matrices
• To perform manipulations such as multiplication or inversion, first convert the matrix to either CSC or CSR format.
• All conversions among the CSR, CSC, and COO formats are efficient, linear-time operations.

Today: CT scans¶

"Can Maths really save your life? Of course it can!!" (lovely article)

(CAT and CT scan refer to the same procedure. CT scan is the more modern term)

This lesson is based off the Scikit-Learn example Compressive sensing: tomography reconstruction with L1 prior (Lasso)

Our goal today¶

Take the readings from a CT scan and construct what the original looks like.

For each x-ray (at a particular position and particular angle), we get a single measurement. We need to construct the original picture just from these measurements. Also, we don't want the patient to experience a ton of radiation, so we are gathering less data than the area of the picture.

Review¶

In the previous lesson, we used Robust PCA for background removal of a surveillance video. We saw that this could be written as the optimization problem:

$$minimize\; \lVert L \rVert_* + \lambda\lVert S \rVert_1 \\ subject\;to\; L + S = M$$

Question: Do you remember what is special about the L1 norm?

Today¶

We will see that:

Imports¶

In [1]:
%matplotlib inline
import numpy as np, matplotlib.pyplot as plt, math
from scipy import ndimage, sparse

In [2]:
np.set_printoptions(suppress=True)


Generate Data¶

Intro¶

We will use generated data today (not real CT scans). There is some interesting numpy and linear algebra involved in generating the data, and we will return to that later.

Code is from this Scikit-Learn example Compressive sensing: tomography reconstruction with L1 prior (Lasso)

Generate pictures¶

In [3]:
def generate_synthetic_data():
rs = np.random.RandomState(0)
n_pts = 36
x, y = np.ogrid[0:l, 0:l]
mask_outer = (x - l / 2) ** 2 + (y - l / 2) ** 2 < (l / 2) ** 2
mx,my = rs.randint(0, l, (2,n_pts))
return res ^ ndimage.binary_erosion(res)

In [208]:
l = 128
data = generate_synthetic_data()

In [209]:
plt.figure(figsize=(5,5))
plt.imshow(data, cmap=plt.cm.gray);


What generate_synthetic_data() is doing¶

In [155]:
l=8; n_pts=5
rs = np.random.RandomState(0)

In [156]:
x, y = np.ogrid[0:l, 0:l]; x,y

Out[156]:
(array([[0],
[1],
[2],
[3],
[4],
[5],
[6],
[7]]), array([[0, 1, 2, 3, 4, 5, 6, 7]]))
In [170]:
x + y

Out[170]:
array([[ 0,  1,  2,  3,  4,  5,  6,  7],
[ 1,  2,  3,  4,  5,  6,  7,  8],
[ 2,  3,  4,  5,  6,  7,  8,  9],
[ 3,  4,  5,  6,  7,  8,  9, 10],
[ 4,  5,  6,  7,  8,  9, 10, 11],
[ 5,  6,  7,  8,  9, 10, 11, 12],
[ 6,  7,  8,  9, 10, 11, 12, 13],
[ 7,  8,  9, 10, 11, 12, 13, 14]])
In [157]:
(x - l/2) ** 2

Out[157]:
array([[ 16.],
[  9.],
[  4.],
[  1.],
[  0.],
[  1.],
[  4.],
[  9.]])
In [59]:
(x - l/2) ** 2 + (y - l/2) ** 2

Out[59]:
array([[ 32.,  25.,  20.,  17.,  16.,  17.,  20.,  25.],
[ 25.,  18.,  13.,  10.,   9.,  10.,  13.,  18.],
[ 20.,  13.,   8.,   5.,   4.,   5.,   8.,  13.],
[ 17.,  10.,   5.,   2.,   1.,   2.,   5.,  10.],
[ 16.,   9.,   4.,   1.,   0.,   1.,   4.,   9.],
[ 17.,  10.,   5.,   2.,   1.,   2.,   5.,  10.],
[ 20.,  13.,   8.,   5.,   4.,   5.,   8.,  13.],
[ 25.,  18.,  13.,  10.,   9.,  10.,  13.,  18.]])
In [60]:
mask_outer = (x - l/2) ** 2 + (y - l/2) ** 2 < (l/2) ** 2; mask_outer

Out[60]:
array([[False, False, False, False, False, False, False, False],
[False, False,  True,  True,  True,  True,  True, False],
[False,  True,  True,  True,  True,  True,  True,  True],
[False,  True,  True,  True,  True,  True,  True,  True],
[False,  True,  True,  True,  True,  True,  True,  True],
[False,  True,  True,  True,  True,  True,  True,  True],
[False,  True,  True,  True,  True,  True,  True,  True],
[False, False,  True,  True,  True,  True,  True, False]], dtype=bool)
In [61]:
plt.imshow(mask_outer, cmap='gray')

Out[61]:
<matplotlib.image.AxesImage at 0x7efcd9303278>
In [62]:
mask = np.zeros((l, l))
mx,my = rs.randint(0, l, (2,n_pts))

Out[62]:
array([[ 0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
[ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.],
[ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.],
[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
[ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.]])
In [63]:
plt.imshow(mask, cmap='gray')

Out[63]:
<matplotlib.image.AxesImage at 0x7efcd9293940>
In [64]:
mask = ndimage.gaussian_filter(mask, sigma=l / n_pts)

In [65]:
plt.imshow(mask, cmap='gray')

Out[65]:
<matplotlib.image.AxesImage at 0x7efcd922c0b8>
In [66]:
res = np.logical_and(mask > mask.mean(), mask_outer)
plt.imshow(res, cmap='gray');

In [67]:
plt.imshow(ndimage.binary_erosion(res), cmap='gray');

In [68]:
plt.imshow(res ^ ndimage.binary_erosion(res), cmap='gray');


Generate Projections¶

Code¶

In [72]:
def _weights(x, dx=1, orig=0):
x = np.ravel(x)
floor_x = np.floor((x - orig) / dx)
alpha = (x - orig - floor_x * dx) / dx
return np.hstack((floor_x, floor_x + 1)), np.hstack((1 - alpha, alpha))

def _generate_center_coordinates(l_x):
X, Y = np.mgrid[:l_x, :l_x].astype(np.float64)
center = l_x / 2.
X += 0.5 - center
Y += 0.5 - center
return X, Y

In [73]:
def build_projection_operator(l_x, n_dir):
X, Y = _generate_center_coordinates(l_x)
angles = np.linspace(0, np.pi, n_dir, endpoint=False)
data_inds, weights, camera_inds = [], [], []
data_unravel_indices = np.arange(l_x ** 2)
data_unravel_indices = np.hstack((data_unravel_indices,
data_unravel_indices))
for i, angle in enumerate(angles):
Xrot = np.cos(angle) * X - np.sin(angle) * Y
inds, w = _weights(Xrot, dx=1, orig=X.min())
mask = (inds >= 0) & (inds < l_x)
camera_inds += list(inds[mask] + i * l_x)
proj_operator = sparse.coo_matrix((weights, (camera_inds, data_inds)))
return proj_operator


Projection operator¶

In [210]:
l = 128

In [211]:
proj_operator = build_projection_operator(l, l // 7)

In [212]:
proj_operator

Out[212]:
<2304x16384 sparse matrix of type '<class 'numpy.float64'>'
with 555378 stored elements in COOrdinate format>

dimensions: angles (l//7), positions (l), image for each (l x l)

In [213]:
proj_t = np.reshape(proj_operator.todense().A, (l//7,l,l,l))


The first coordinate refers to the angle of the line, and the second coordinate refers to the location of the line.

The lines for the angle indexed with 3:

In [214]:
plt.imshow(proj_t[3,0], cmap='gray');

In [215]:
plt.imshow(proj_t[3,1], cmap='gray');

In [216]:
plt.imshow(proj_t[3,2], cmap='gray');

In [217]:
plt.imshow(proj_t[3,40], cmap='gray');


Other lines at vertical location 40:

In [218]:
plt.imshow(proj_t[4,40], cmap='gray');

In [219]:
plt.imshow(proj_t[15,40], cmap='gray');

In [220]:
plt.imshow(proj_t[17,40], cmap='gray');


Intersection between x-rays and data¶

Next, we want to see how the line intersects with our data. Remember, this is what the data looks like:

In [221]:
plt.figure(figsize=(5,5))
plt.imshow(data, cmap=plt.cm.gray)
plt.axis('off')
plt.savefig("images/data.png")

In [222]:
proj = proj_operator @ data.ravel()[:, np.newaxis]


An x-ray at angle 17, location 40 passing through the data:

In [223]:
plt.figure(figsize=(5,5))
plt.imshow(data + proj_t[17,40], cmap=plt.cm.gray)
plt.axis('off')
plt.savefig("images/data_xray.png")


Where they intersect:

In [224]:
both = data + proj_t[17,40]
plt.imshow((both > 1.1).astype(int), cmap=plt.cm.gray);


The intensity of an x-ray at angle 17, location 40 passing through the data:

In [225]:
np.resize(proj, (l//7,l))[17,40]

Out[225]:
6.4384498372605989

The intensity of an x-ray at angle 3, location 14 passing through the data:

In [226]:
plt.imshow(data + proj_t[3,14], cmap=plt.cm.gray);


Where they intersect:

In [227]:
both = data + proj_t[3,14]
plt.imshow((both > 1.1).astype(int), cmap=plt.cm.gray);


The measurement from the CT scan would be a small number here:

In [228]:
np.resize(proj, (l//7,l))[3,14]

Out[228]:
2.1374953737965541
In [229]:
proj += 0.15 * np.random.randn(*proj.shape)


In [230]:
a = [1,2,3]
b= [4,5,6]

In [231]:
c = list(zip(a, b))

In [232]:
c

Out[232]:
[(1, 4), (2, 5), (3, 6)]
In [233]:
list(zip(*c))

Out[233]:
[(1, 2, 3), (4, 5, 6)]

In [234]:
plt.figure(figsize=(7,7))
plt.imshow(np.resize(proj, (l//7,l)), cmap='gray')
plt.axis('off')
plt.savefig("images/proj.png")


Regresssion¶

Now we will try to recover the data just from the projections (the measurements of the CT scan)

Linear Regression: $Ax = b$¶

Our matrix $A$ is the projection operator. This was our 4d matrix above (angle, location, x, y) of the different x-rays:

In [203]:
plt.figure(figsize=(12,12))
plt.title("A: Projection Operator")
plt.imshow(proj_operator.todense().A, cmap='gray')

Out[203]:
<matplotlib.image.AxesImage at 0x7efcd4199a20>

We are solving for $x$, the original data. We (un)ravel the 2D data into a single column.

In [202]:
plt.figure(figsize=(5,5))
plt.title("x: Image")
plt.imshow(data, cmap='gray')

plt.figure(figsize=(4,12))
# I am tiling the column so that it's easier to see
plt.imshow(np.tile(data.ravel(), (80,1)).T, cmap='gray')

Out[202]:
<matplotlib.image.AxesImage at 0x7efcd3903be0>

Our vector $b$ is the (un)raveled matrix of measurements:

In [238]:
plt.figure(figsize=(8,8))
plt.imshow(np.resize(proj, (l//7,l)), cmap='gray')

plt.figure(figsize=(10,10))
plt.imshow(np.tile(proj.ravel(), (20,1)).T, cmap='gray')

Out[238]:
<matplotlib.image.AxesImage at 0x7efcd3bebf28>

Scikit Learn Linear Regression¶

In [84]:
from sklearn.linear_model import Lasso
from sklearn.linear_model import Ridge

In [126]:
# Reconstruction with L2 (Ridge) penalization
rgr_ridge = Ridge(alpha=0.2)
rgr_ridge.fit(proj_operator, proj.ravel())
rec_l2 = rgr_ridge.coef_.reshape(l, l)
plt.imshow(rec_l2, cmap='gray')

Out[126]:
<matplotlib.image.AxesImage at 0x7efcd453d5c0>
In [179]:
18*128

Out[179]:
2304
In [ ]:
18 x 128 x 128 x 128

In [178]:
proj_operator.shape

Out[178]:
(2304, 16384)
In [87]:
# Reconstruction with L1 (Lasso) penalization
# the best value of alpha was determined using cross validation
# with LassoCV
rgr_lasso = Lasso(alpha=0.001)
rgr_lasso.fit(proj_operator, proj.ravel())
rec_l1 = rgr_lasso.coef_.reshape(l, l)
plt.imshow(rec_l1, cmap='gray')

Out[87]:
<matplotlib.image.AxesImage at 0x7efcd4919cf8>

The L1 penalty works significantly better than the L2 penalty here!