DCAN Lab Introduction to Image Mathematics¶

This is written as introductory material to the mathematics of images for neuroscientists and psychologists interested in neuroimaging. It covers some basic terminology, image arithmetic, as well as some more advanced topics relevant to medical imaging.

Topics¶

• What makes an image?

• What are some operations on images?

  • Images and Numbers

  • Images and other Images

  • Masking and Segmentation

  • Transformations

• What kind of images are there in MRI research?

  • Anatomical

  • Functional

  • Field Maps

  • DWI

• Common MRI Processing Tools

  • fsl

  • FreeSurfer

  • ANTs

  • Workbench Command

Interactive Setup¶

If you wish to use this notebook interactively, there are a number of python packages you will need to install first. The very fist cell contains all of the import statements necessary. Install anaconda or setup a virtual environment, and ensure that all packages in the requirements.txt file are installed.

What Makes An Image?¶

In this context we will always be talking about digital images, the way that they are represented on a computer. Starting with 2-dimensional images, a standard digital image is represented on a computer by a grid or matrix of numbers.

If you are new to Jupyter notebooks such as this one, the off-white boxes such as the one below, with In [N]: (where N is some number) in the margins, are the python code used to display images/examples. The code may be followed or ignored for the purposes of this introduction, depending on your interest in python/programming.

Above you can see the grids we have created, which we will call grid1 and grid2. Notice grid1 is split half and half between 1 and 0, grid2 is alternating 1s and 0s.

Above are text representations of numpy arrays, they basically follow the shape of their corresponding matrices, which we'll see below.

What you see is organized by two indices: a row index and a column index which together map each element.

Typically, we use the index variables i and j as subscripts to correspond to rows and columns respectively.

For a matrix, $A$, we use $a_{ij}$ to represent each individual value in the matrix (commas are generally omitted in single-digit cases):

$$\begin{bmatrix}a_{1,1} & a_{1,2} & a_{1,3} & a_{1,4} & a_{1,5} & a_{1,6}\\a_{2,1} & a_{2,2} & a_{2,3} & a_{2,4} & a_{2,5} & a_{2,6}\\a_{3,1} & a_{3,2} & a_{3,3} & a_{3,4} & a_{3,5} & a_{3,6}\\a_{4,1} & a_{4,2} & a_{4,3} & a_{4,4} & a_{4,5} & a_{4,6}\\a_{5,1} & a_{5,2} & a_{5,3} & a_{5,4} & a_{5,5} & a_{5,6}\\a_{6,1} & a_{6,2} & a_{6,3} & a_{6,4} & a_{6,5} & a_{6,6}\end{bmatrix}$$

In grid above, the $a_{ij}$ variables corresponds to the real valued numbers in the grids below:

$$\begin{bmatrix}1 & 1 & 1 & 0 & 0 & 0\\1 & 1 & 1 & 0 & 0 & 0\\1 & 1 & 1 & 0 & 0 & 0\\1 & 1 & 1 & 0 & 0 & 0\\1 & 1 & 1 & 0 & 0 & 0\\1 & 1 & 1 & 0 & 0 & 0\end{bmatrix} \quad \begin{bmatrix}1 & 0 & 1 & 0 & 1 & 0\\0 & 1 & 0 & 1 & 0 & 1\\1 & 0 & 1 & 0 & 1 & 0\\0 & 1 & 0 & 1 & 0 & 1\\1 & 0 & 1 & 0 & 1 & 0\\0 & 1 & 0 & 1 & 0 & 1\end{bmatrix}$$

Both images and matrices are grids of points (or pixels), so we can use this same notation for images.

Let's see how these matrices look as a images:

Notice that wherever we had a 1, the image shows a white square. Wherever we had a 0, the image shows a black square. Whether the number or the visible representation, these are pixels.

If we were discussing photographs such as taken on a digital camera, each pixel would actually have three values associated: R, G, and B. In the context of MRI, we are usually only concerned with grayscale images (colloquially, black and white images).

The images above, which only have two values: 0, and 1, are a special case of greyscale image, known as a binary image

What Are Some Operations on Images?¶

Here we will discuss mathematical operations which can be performed on images.

Images and Numbers¶

Given some number $x$, and another number $a$, I can write some simple equations like:

$$ x + a = 5 \\ x - a = 1 \\ x \times a = 6 \\ x ^ a = 9 $$

Which you will recognize as addition, subtraction, multiplication and exponentiation.

For some image $I$, and a number $a$, I can carry out all of the same operations, for example, I can add $a$ to grid1 (the half-and-half grid from before):

$$ I + a = ?$$

As it turns out, the extension is quite simple.

Notice that the value of each element in grid1 has increased by 2! We call these elementwise operations, where the normal mathematical operation: addition, multiplication, etc. is carried out individually on each pixel/element of the image:

$$ J_{n, m} = a + I_{n, m} $$

where $I_{n, m}$ is the element in row n, column m, of $I$

Adding 2 to an image is the same as adding 2 to every pixel, multiplying an image by 2 is the same as multiplying every pixel by 2. In fact, even fancy operations like log, exponents, or trig functions obey this same rule.

Let's take a look at the visual result for illustrative purposes:

Question: Why does the image still look the same as before?¶

...

Answer: You must be careful when visually comparing different images to make sure they are scaled the same. 2 is here black because it is the lowest value in the image. 3 is white because it is the highest. We could also scale this image from 0 to 4, and we would see some slightly different shades of gray in place of 3 and 2:

Images and other Images¶

Similar to before, we can also have addition, subtraction, multiplication, division, etc. between images. For example, given an image $I$ and an image $J$ we can still have:

$I + J = K$

where K is a new image. To perform this operation, we add I and J elementwise, like so:

$$ K_{n,m} = I_{n,m} + J_{n,m} $$

Here we will add the two images from above, the checkered and the half and half:

Note that there are now three distinct values in the matrix: 2, 1, 0, and three distinct gray levels in the resulting image on the right.

Try it for yourself!¶

Code cells can be run by hitting [Shift] + [ENTER]

I'm going to modify the code a bit here, the visual minimum and maximum will be whatever the min and max of I, J, and K are together.

Some Real Images¶

okay, so that's all well and good but these images are kind of boring, let's spice it up.

Let's see what happens when we multiply these two together!

In the cell below, write the equation for multiplying the images in place of the None, then hit [Shift] + [Enter] to run the code.

Note: For interactive use of this notebook, you will need to set up a proper environment with all of the required packages.

Transformations¶

Now we will start to get into the idea of transforming an image. Transformations are what allow us to register images.

Linear Transformations¶

There are 4 aspects of linear transformation:

• Translation

• Rotation

• Scaling

• Shearing

Affine Linear Transforms¶

Below is an illustration of each of the 4, sometimes known by slightly different names. Here "Shearing" is "Skew", and "Differential" Scaling is used because the scaling need not be equivilent in different dimentions (hence the differential), turning the square into a rectangle.

Together these 4 transformation types make up "affine" transformations. A single linear transformation may involve one, or all of them. It is specifically called 'affine' when we are using or allowed to use all 4.

Affine transformations are the set of continuous transformations which preserve parallel lines. It is helpful to visualize a square when thinking about these transforms, because it beautifully illustrates the different possibilities: Any transformation of a square into a rectangle or rhombus is permitted through skew and scaling, but trapezoids and other quadrilaterals would not be possible, as it would break the parallelism.

When you hear the term degrees of freedom used in the context of affine transformations, it corresponds to the parameters which dictate these 4 types of transforms.

By constraining certain parameters, we can get a special classes of linear transformation:

Rigid Body Linear Transforms¶

Translation and rotation, when considered alone, encompass rigid body or euclidean transformations. A translation is a shift of the object from point A to point B. A rotation is performed about the origin. Transformations involving only these two aspects are called "rigid body", because the shape and size of the object are constant through these transforms.

Below you can use the sliders to play with the translation and rotation. Notice what happens to the edges of the square as it moves through discrete pixels in the image's field of view.

You may notice a couple of things here:

1. the translation tx and ty are inverted

   this is akin to shifting algebraic functions f(x) to the left or right on a graph.

1. the origin is in the top-left instead of the bottom-left

   this is how image coordinates differ from cartesian coordinates.

1. the square rotates about the origin, not its own center of mass.

   This is because the transformation is on the entire coordinate system. If there were two objects, they would both be transformed according to the coordinate transform.

In order to flip an object upside down, you need to rotate it 180 degrees about the origin, then translate it back across (twice its starting distance from $(0,0)$. Later, when we get into MRI images we will talk about sform and qform, and why they help with this process.

You should also note that the edges blur as the square moves around in space. This is because the edges of the square do not land neatly into a single pixel.

Interpolation¶

Bi-Linear¶

The default option on this interactive widget, where the blurring is observed, is called "Bi-Linear" interpolation. New pixels are filled in with the average value of the pixels around them. That way we get a sense of where the square has 'mostly' entered a pixel, or 'slightly' entered a pixel by the resulting darker fade.

Nearest Neighbor¶

One option would be a sort of "majority vote": any pixel which is closer to a square edge than empty space gets filled in. Click the "Nearest Neighbor" box to see this in action. You can see how the shape can lose its squarish form with this method. To begin with, the square is intensity 1 and the background is 0. For the faded pixels from before, any with a value less than .5 become 0 and any greater than .5 become 1.

Q: Why use Nearest Neighbor if the result can be so mishapen?

A: It is primarily useful when dealing with labeled data. If the square/object is a box dileaneating someone's hippocampus, faded values may not have any discernable meaning. Furthermore, the faded values less than .5 are actually mostly "not square", so depending on how specific we want to be, it may be important to not include those points.

Interactive Affine Transformations¶

Below, I will allow you to play with affine transformations. You now have an additional 4 parameters to play with while moving the image. See what they do!

(hint: $k_x$ and $k_y$ may appear to do nothing at first, but try playing with $\theta$ as well!)

If the square gets lost, press the reset button.

3 dimensional images (MRI)¶

without further ado, let's take a look at these same concepts on 3 dimensional brains! The concepts of affine transformations generalizes nicely here. In fact, you probably have more of an intuitive grasp of them, as your mind is perfectly suited to modeling 3D transformations.

The only trouble is that it can be somewhat tedious to constantly visualize what the computer is doing. A lot of neuroscience researchers use a lot of transformation tools such as fsl flirt without a good intuition for the rules of what its doing. By playing with some of it interactively, we might get a better intuition for how they work, and what they're good for.

Now, remember rigid body transformations? Let's revisit them on an MRI image. Bear with me for a second, this will take a little bit of code to setup...

Notice that we now have 3 angles about which we can rotate, instead of the single rotation axis from before. They are called "euler angles", and typically represented by greek letters, but can also be referred to in aircraft terms, "yaw", "pitch", and "roll".

This is also a great opportunity to revisit interpolation. I have provided four available options:

• Nearest Neighbor
• Linear
  • for 2d images, it is also called "bilinear interpolation"
  • for 3d images, "trilinear"
• Quadratic
  • also called "thin-plate spline interpolation"
• Cubic
  • also called "cubic spline interpolation" or sometimes just "spline interpolation"

Try rotating/translating the image, and then changing the interpolation order.

See if you can answer these questions:

1. why does interpolation change nothing if all of the parameters are set to zero?

1. now turn a few dials. Try to find a place where "Nearest Neighbor" looks especially blocky. Now play with the interpolation. What do you observe? Which method would you use here?

1. now reset to zero, and rotate the image by 90 degrees for any angle. Try changing the interpolation. What is happening?

To answer the last question explicitly, if a transform is such that the "new" pixels of the image fall exactly on top of the old pixels, no interpolation is needed! In practice, for computed registrations, this will seldom if ever be the case.

So, for 90 degree rotations, and exact translations by integer numbers of pixels, there is no interpolation needed. This is relevant for "reorientation" of images, where we can rest assured that the data remains untouched.

Nonlinear Transformations¶

To be continued...