Understanding the maths of Computed Tomography (CT) scans

January 9, 2019
By

(This article was first published on R-Bloggers – Learning Machines, and kindly contributed to R-bloggers)


Noseman is having a headache and as an old-school hypochondriac he goes to see his doctor. His doctor is quite worried and makes an appointment with a radiologist for Noseman to get a CT scan.

Modern CT scanner from Siemens

Because Noseman always wants to know how things work he asks the radiologist about the inner workings of a CT scanner.

The basic idea is that X-rays are fired from one side of the scanner to the other. Because different sorts of tissue (like bones, brain cells, cartilage etc.) block different amounts of the X-rays the intensity measured on the other side varies accordingly.

The problem is of course that a single picture cannot give the full details of what is inside the body because it is a combination of different sorts of tissue in the way of the respective X-rays. The solution is to rotate the scanner and combine the different slices.

How, you ask? Good old linear algebra to the rescue!

We start with the initial position and fire X-rays with an intensity of 30 (just a number for illustrative purposes) through the body:

Initial position

As can be seen in the picture the upper ray goes through areas 1, 2 and 3 and let’s say that the intensity value of 12 is measured on the other side of the scanner:

    \[30-x_1-x_2-x_3=12\]

or

    \[x_1+x_2+x_3=18\]

The rest of the formula is found accordingly:

    \[\underbrace{ \begin{pmatrix}   1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\   0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 \\   0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 \end{pmatrix} }_{\bold{A}_1} \underbrace{ \begin{pmatrix}   x_1 \\   x_2 \\   x_3 \\   x_4 \\   x_5 \\   x_6 \\   x_7 \\   x_8 \\   x_9 \end{pmatrix} }_{\bold{x}} = \underbrace{ \begin{pmatrix}   18 \\   21 \\   18 \end{pmatrix} }_{\bold{b}_1}\]

We then rotate the scanner for the first time…

Position after first rotation

…which gives the following formula:

    \[\underbrace{ \begin{pmatrix}   0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 \\   0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 \\   0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 \end{pmatrix} }_{\bold{A}_2} \underbrace{ \begin{pmatrix}   x_1 \\   x_2 \\   x_3 \\   x_4 \\   x_5 \\   x_6 \\   x_7 \\   x_8 \\   x_9 \end{pmatrix} }_{\bold{x}} = \underbrace{ \begin{pmatrix}   18 \\   21 \\    9 \end{pmatrix} }_{\bold{b}_2}\]

And a second rotation…

Position after second rotation

…yields the following formula:

    \[\underbrace{ \begin{pmatrix}   0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\   0 & 0 & 1 & 0 & 1 & 0 & 1 & 0 & 0 \\   0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 \end{pmatrix} }_{\bold{A}_3} \underbrace{ \begin{pmatrix}   x_1 \\   x_2 \\   x_3 \\   x_4 \\   x_5 \\   x_6 \\   x_7 \\   x_8 \\   x_9 \end{pmatrix} }_{\bold{x}} = \underbrace{ \begin{pmatrix}   18 \\   14 \\   16 \end{pmatrix} }_{\bold{b}_3}\]

Now we are combining all three systems of equations:

    \[\begin{pmatrix}   \bold{A}_1 \\   \bold{A}_2 \\   \bold{A}_3 \end{pmatrix} \bold{x} = \begin{pmatrix}   \bold{b}_1 \\   \bold{b}_2 \\   \bold{b}_3 \end{pmatrix}\]

or written out in full:

    \[\underbrace{ \begin{pmatrix}   1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\   0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 \\   0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 \\   0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\   0 & 0 & 1 & 0 & 1 & 0 & 1 & 0 & 0 \\   0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 \\   0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\   0 & 0 & 1 & 0 & 1 & 0 & 1 & 0 & 0 \\   0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 \end{pmatrix} }_{\bold{A}} \underbrace{ \begin{pmatrix}   x_1 \\   x_2 \\   x_3 \\   x_4 \\   x_5 \\   x_6 \\   x_7 \\   x_8 \\   x_9 \end{pmatrix} }_{\bold{x}} = \underbrace{ \begin{pmatrix}   18 \\   21 \\   18 \\   18 \\   21 \\    9 \\   18 \\   14 \\   16 \end{pmatrix} }_{\bold{b}}\]

Here is the data of the matrix \bold{A} for you to download: ct-scan.txt).

We now have 9 equations with 9 unknown variables… which should easily be solvable by R, which can also depict the solution as a gray-scaled image… the actual CT-scan!

A <- read.csv("data/ct-scan.txt")
b <- c(18, 21, 18, 18, 21, 9, 18, 14, 16)
v <- solve(A, b)
matrix(v, ncol = 3, byrow = TRUE)
##      [,1] [,2] [,3]
## [1,]    9    9    0
## [2,]    9    5    7
## [3,]    9    9    0
image(matrix(v, ncol = 3), col = gray(4:0 / 4))
CT of Noseman

The radiologist looks at the picture… and has good news for Noseman: everything is how it should be! Noseman is relieved and his headache is much better now…

Real CT scans make use of the same basic principles (of course with a lot of additional engineering and maths magic 😉 )

Here are real images of CT scans of a human brain…

Source: Wikimedia

… which can be combined into a 3D-animation:

Source: Wikimedia

Isn’t it fascinating how a little bit of maths can save lives!

To leave a comment for the author, please follow the link and comment on their blog: R-Bloggers – Learning Machines.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Search R-bloggers


Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)