Image processing is a very useful tool for scientists in the lab, and for everyday uses as well. For, example, an astronomer may use image processing to help find and recognize stars, or a self driving car may use it to stay in the correct lane.
This lecture will teach you some of the fundamental tools of image processing so that you can use it in your own studies/applications.
Figure 1: Felix the cat
As we can see from Felix the cat, images are just a coordinate system in which each block (x,y) has a value associated to it. These values are later interpreted as a color.
For example, in Figure 1, we would be describing two colors (black and white) and we could describe them using 1's and 0's. In general, 0 is taken to mean the absence of color, which means that 0 = black and 1 = white. Putting these ideas together, we can see that each point in an image requires three peices of information:
So if we were instructing a computer to produce even a small image, we would need to give it a large list of numbers.
Figure 2: Checker board pattern
For example, suppose we wanted the computer to produce the checkboard patter that we see in Figure 2. If we make a list of the information the computer needs, we can format it like (x-coordinate, y-coordinate, color) and the whole list would be
(1,1,1), (1,2,0), (1,3,1), (2,1,0), (2,2,1), (2,3,0), (3,1,1), (3,2,0), (3,3,1)
Figure 3: Images
To make the image making process easier, people decided to ditch the traditional coordinate system and use matricies instead! However, because both systems work and because sometimes one method can be more convenient than the other, both still exists, but they have different names.
When you use a coordinate system to instruct the computer, that's a scatter plot. When you use a matrix, that's called an image.
To instruct the computer to make a scatter plot of the 9 square checkerboard, we had to give it 3x9=27 numbers. To insturct it to make an image of the checkerboard, we only have to give it 9 numbers, but they have to be in the correct order. Specifically, the order looks like this:
image = [ [1, 0, 1],
[0, 1, 0],
[1, 0, 1] ]Each location that we assign a value to is considered a pixel. Thus, we have created an image that is 3 pixels wide and 3 pixels tall, and it has a total of 9 pixels.
When we watch youtube videos at 1080p, we are actually looking at pictures that are 1080 pixels tall and 1920 pixel wide. These images have a total of 1080x1920 = 2,073,600 pixels. If we round to the nearest million, there are approximately 2 million pixels in each image. This would be considered a 2 Mega Pixel (MP) image. This is the same number phone manufactures using when talking about how many megapixels their newest device has.
In [1]:
# set exercise1 equal to your matrix
exercise1 = #your matrix goes here
We can check if we got this right by having the computer make the image. To do this, we will use matplotlib's "plt.imshow" function, which takes in matricies and turns them into images.
In [2]:
import matplotlib.pyplot as plt
import numpy as np
plt.imshow(exercise1)
Out[2]:
Now, we have gotten a figure that looks like the correct board, but the colors are wrong. To fix this, we need to tell the computer that we want to see the image in grayscale. We do this by inserting a "cmap" into the plt.imshow function. "cmap" stands for "color map" and in this case, we will set it to grayscale
In [3]:
from matplotlib import cm
plt.imshow(exercise1, cmap=cm.gray)
Out[3]:
If the output of the code above does not match the chekerboad in the exercise try again or ask for help. Notice how the order of the numbers is important to final image.
Now we will use a scatter plot to create the checkerboard. This time we will use matplotlib's "plt.scatter" function to create the figure.
In the code below, we have given you a list a coordinates and values for the blocks at those coordinates. Correct the list so that it reproduces the checkerboard from Exercise 1.
In [4]:
exercise2 = [[1,1,0], [1,2,1], [1,3,0], [2,1,1], [2,2,0], [2,3,0], [3,1,1], [3,2,0], [3,3,1]]
exercise2 = np.array(exercise2)
x_coordinates = exercise2[::,0]
y_coordinates = exercise2[::,1]
values = exercise2[::,2]
plt.scatter(x_coordinates, y_coordinates, c=values, marker='s',s = 500, cmap=cm.gray)
Out[4]:
Recall that the order of the numbers doesn't matter for scatter plots. To see this, go back to the previous code block, and switch the order of the exercise2 list. Even though you change the order, the image produced doesn't change.
Now that we have covered the fundamentals, we can start diving into image manipulations starting with slices and arithmatic.
There are many tools we can use to learn about image processing. Today, we will be using Scikit-Image, a free to use and easy to learn python library, and we'll start by looking at an image of an insect.
In [6]:
# first we import the input/output module from skimage called "io". io lets us
# read files on the computer or online
from skimage import io
insect = io.imread('https://matplotlib.org/3.1.1/_images/stinkbug.png')
#we can take a look at the image we just imported once again using image show.
plt.imshow(insect, cmap=cm.gray, vmin=0, vmax=256)
plt.colorbar()
Out[6]:
Remember that, since we are working with images, the information is stored inside of a matrix. Lets take a look at that matrix.
In [7]:
# show the array
print(insect)
The first thing we notice is that this image is no longer just 1's and 0's. The numbers go from 0 to 255, and recalling that 0 = absence of light, we take 0=black and 255=white. All of the numbers between 0 and 255 are various shades of gray.
Next, we can see that there are ellipses inside of the matrix. This is python's way of telling us that there are so many positions and values that writing them all out would take a huge amount of space on the screen, and a list that big is hard to read.
In the checker board example we had an image that had 3 rows and 3 columns. We can check the shape of any matrix by adding ".shape" after its name. Let's see the shape of the insect image
In [8]:
# show the shape of the array
print(insect.shape)
It looks like this matrix has 375 rows and 500 columns, so we expect it to be wider than it is tall. We can take a look at the value of a specific pixel by calling its location like this:
In [9]:
# show the value of the pixel in the 100th row and 200th column
print(insect[100,200])
Cropping
Cropping an image is very simple. Basically, we have a really big matrix, but want to get rid of the rows and columns that we are not interestesd in. We can do this by slicing the matrix usng brackets [ ] as show below:
In [10]:
# display only the section of the image between columns 150 - 350 and rows 200 - 400
plt.imshow(insect[50:150,150:275], cmap=cm.gray, vmin=0, vmax=256)
Out[10]:
In this case we cropped to the left antenna.
In [11]:
plt.imshow(insect[insert_your_slice_here], cmap=cm.gray, vmin=0, vmax=256)
Out[11]:
In [12]:
#divide the image by two and display it
plt.imshow(insect/2, cmap=cm.gray, vmax=255, vmin=0)
plt.colorbar()
Out[12]:
To brighten the image, we can instead multiply the matrix by 1.5:
In [13]:
#multiply the image by 1.5 and display it
plt.imshow(insect*1.5, cmap=cm.gray, vmax=255, vmin=0)
plt.colorbar()
Out[13]:
We can also add and subtract images from one another. For example, we can subtract an image from its self:
In [14]:
# subtract the image to itself and display
plt.imshow(insect-insect, cmap=cm.gray, vmax=256, vmin=0)
# check that the array is full of zeros
print(insect-insect)
We can also directly modify the values of the matrix elements ourselves. For example, if we want to set a portion of the image to solid black, we can do so by setting all of the pixels in that region to zero, like this:
In [15]:
# set the region x=[50,100] and y=[50,100] to zero and display
img = np.copy(insect)
img[50:100,50:100] = 0
plt.imshow(img, cmap=cm.gray, vmax=256, vmin=0)
Out[15]:
In [16]:
img = np.copy(insect)
img['insert_slice'] = 'insert value'
plt.imshow(img, cmap=cm.gray, vmax=256, vmin=0)
Out[16]:
In [17]:
# challenge problem: make a fading image from top to bottom (hint: multiply the image by a gradient you can create)
Out[17]:
Advanced Image Processing
In essence, "blob detection" is the process of searching for and identifying bright spots on a dark background, or dark spots on a light background. The "blobs" can correspond to any number of physical phenomena. In astronomy, for example, this can be relevant if you have an image of a part of the sky and want to identify something like galaxies. You might be wondering though, why not just look at an image yourself and pick out blobs "by eye"? In our example of identifying galaxies in an astronomical image, a systematic way of doing this may be beneficial for several reasons:
These are just some of the reasons why a systematic approach to identifying blobs in images could be beneficial. Checking the output of your algorithm by eye, however, is wise to make sure that it is not outputting nonsense.
Let's get started with first reading in our first astronomical image; a nearby group of galaxies! Astronomical images are stored as "fits" files, this is essentially a matrix that contains the information of how bright the sky is at each pixel in your image, and a "header", which contains information about how pixels translate into coordinates on the sky, what instrument was used, and the units of the image for example. We will read the fits file into a numpy array using the package called astropy (see http://docs.astropy.org/en/stable/index.html for documentation). This package streamlines the process of working with astronomical images in Python.
First we will import a few packages that we are going to use throughout this section: skimage will allow us to run the blob detection algorithms, matplotlib will allow us to plot our data, astropy will allow us to read "fits" files, and numpy will let us do cool things with arrays.
In [18]:
from skimage.feature import blob_dog, blob_log, blob_doh
import matplotlib.pyplot as plt
from astropy.io import fits
import numpy as np
This line reads the fits file --> data table and header.
In [19]:
hdu = fits.open('./RSCG1.fits')
In [85]:
#Insert your code here:
The following line reads the data associated with the header into a numpy array. A fits file can have more than one header, and since Python indices start at 0, we tell Python to read the data associated with the first (and only) header associated with this fits file.
In [21]:
image = hdu[0].data
# This line closes the fits file because we don't need it anymore; we've already read the data into an array.
hdu.close()
In [22]:
#Insert your code here:
print(type(image))
print(image)
print(image.shape)
print(np.max(image))
print(np.min(image))
That's it, you've now read your astronomical image into Python and can work with the data and analyze it! Let's visualize the data now to get an idea of what it is we are actually dealing with.
In [23]:
# imshow will map the 2D array (our data), and origin='lower' tells imshow to plot (0,0) in the bottom left corner.
plt.imshow(image,origin='lower')
plt.colorbar()
plt.show()
# EXERCISE: What happens if you remove origin='lower'? How does the mapping of the image change?
We can now run some blob finding algorithms on this data set and see what we find and if it makes sense! For this, we will use the scikit-image package which is an image processing pacakge in Python. This package has three algorithms for detecting blobs (all discussed here: https://en.wikipedia.org/wiki/Blob_detection):
This is great, but what does this actually mean? How do they work? What's the difference between them?
This algorithm starts by smoothing (blurring) the entire image with a two-dimensional Gaussian function:
$$\begin{align} g(x,y,\sigma) = \frac{1}{2\pi\sigma^2} e^\left(-\frac{x^2+y^2}{2\sigma^2}\right) \end{align}$$where $(x,y)$ correspond to pixel coordinates on your image, and $\sigma$ represents the "size" or "width" of the Gaussian function. The algorithm then performs a mathematical computation called taking the "Laplacian", however we will not go into the details of that. What this effectively gives you is strong responses for blobs of size $\sqrt{2}\sigma$. This means that the algorithm is most sensitive to detecting blobs that have a similar size to that of the Gaussian function that you smooth your image with in the first place. In order to detect blobs of varying size, you need to use an approach that allows you to smooth your image with Gaussians of varying sizes to try to match the scale of potential blobs. This approach is the most accurate but also the slowest approach, especially for the largest blobs, because it needs to smooth the image with a larger Gaussian.
The scikit-image function "blob_log" (https://scikit-image.org/docs/dev/api/skimage.feature.html#skimage.feature.blob_log) has this capability. To run it, we simply need to call the function and provide it with the correct input parameters. In the following piece of code, this is what we are going to provide the function with:
The call to this function will output a list of x and y coordinates corresponding to the center of each blob, as well as the "size" of the Gaussian that it used to find that blob. Let's try it out!
In [24]:
# This line calls the blob detection (Laplacian of Gaussian function) on our galaxies image with the parameters that
# we provide
blobs_log = blob_log(image, min_sigma = 2, max_sigma=9, num_sigma=8, threshold=.05)
In [25]:
input("What kind of output is blobs_log?")
input("What are its dimension?")
input("What does the length of blobs_log physically correspond to?")
Out[25]:
In [26]:
# In this line we use the third column to calculate the radius of the blob itself.
blobs_log[:, 2] = blobs_log[:, 2] * np.sqrt(2)
Now we'll plot the blobs onto our image! For this, we will create circles with plt.Circle using the corrdinate and radius information we have from blobs_log. Then we will use add_patch to add that circle to our image. Since the plotting software needs information about the axes of the image in order to know where each circle goes, we use plt.gca() which means "Get Current Axis" and store the axis information into "ax". Then when we add the circle patches to the image, they will appear in the correct places. Finally, since there is more than one blob, we will create the circular patches and add them to the image one-by-one using a while loop.
In [27]:
ax = plt.gca()
ax.imshow(image,origin='lower')
cnt = 0
while cnt < len(blobs_log):
c = plt.Circle((blobs_log[cnt][1], blobs_log[cnt][0]), blobs_log[cnt][2], color='white', linewidth=2, fill=False)
ax.add_patch(c)
cnt = cnt + 1
plt.show()
In [28]:
blobs_log = blob_log(image, min_sigma = 2, max_sigma=9, num_sigma=8, threshold=.05)
### Play with the paramters above this line ###
blobs_log[:, 2] = blobs_log[:, 2] * np.sqrt(2)
ax = plt.gca()
ax.imshow(image,origin='lower')
cnt = 0
while cnt < len(blobs_log):
c = plt.Circle((blobs_log[cnt][1], blobs_log[cnt][0]), blobs_log[cnt][2], color='white', linewidth=2, fill=False)
ax.add_patch(c)
cnt = cnt + 1
plt.show()
Write down your findings below, such as what parameters give you no blobs, how does changing the parameters affect the results, and anything else that you find interesting!
In [29]:
# RECORD YOUR FINDINGS HERE
This algorithm is an approximation of the Laplacian of Gaussian, where rather than actually computing the "Laplacian" of the Gaussian that we mentioned above, DoG approximates this computation by taking the difference of two successively smoothed (blurred) images. Blobs are then detected as bright-on-dark spots in these difference images. This algorithm has the same disadvantage for larger blobs as the Laplacian of Gaussian algorithm.
To run this algorithm, we use the function "blob_dog" (https://scikit-image.org/docs/stable/api/skimage.feature.html#skimage.feature.blob_dog), which takes the same parameters as blob_log, except for num_sigma, because the algorithm needs to take the difference between successive smoothings so it figures out the number of times it needs to smooth the image. This function returns the same type of array as "blob_log" : (y-coord, x-coord, size).
In [86]:
blobs_dog = blob_dog('fill here')
blobs_dog[:, 2] = blobs_dog[:, 2] * np.sqrt(2)
ax = plt.gca()
ax.imshow(image,origin='lower')
cnt = 0
while cnt < len(blobs_dog):
c = plt.Circle((blobs_dog[cnt][1], blobs_dog[cnt][0]), blobs_dog[cnt][2], color='white', linewidth=2, fill=False)
ax.add_patch(c)
cnt = cnt + 1
plt.show()
The final option that scikit-image provides is the Determinant of Hessian algorithm. This algorithm performs a different kind of mathematical operation on the smoothed images, it calculates the Hessian matrix of the image (but we won't go into that), and searches for maxima in this matrix. It it the fastest approach, but it is also less accurate and has trouble finding very small blobs.
To try it out, we will run the function "blob_doh" (https://scikit-image.org/docs/stable/api/skimage.feature.html#skimage.feature.blob_doh), which takes as input the same parameters as the "blob_log" function. Let's try it out!
In [31]:
blobs_doh = blob_doh('fill here')
ax = plt.gca()
ax.imshow(image,origin='lower')
cnt = 0
while cnt < len(blobs_doh):
c = plt.Circle((blobs_doh[cnt][1], blobs_doh[cnt][0]), blobs_doh[cnt][2], color='white', linewidth=2, fill=False)
ax.add_patch(c)
cnt = cnt + 1
plt.show()
In [ ]:
# RECORD YOUR THOUGHTS HERE
Particle Tracking
Now that we have learned several image processing techniques, we are ready to put them together in an application. In physics, there is an experiment called the Monopole Ion Trap that uses electric fields to trap particles as show below.
The basic priciple is that the particle is attracted to and repeled by the rod at the top of the image. In the image above, we have a very stable, well behaved, trap in which the ion just goes back and forth between two locations. However, under certain curcumstances, the particle may start to behave erratically, as shown below.
We would like to track the particle so that we can study its position and velocity.
Lets start by importing the movie of the stable particle into python. We have saved each movie, image by image, inside of folders on the computers. So we will import the stable movie as follows:
In [32]:
# we import the operating system library OS to help us import many files at once
import os
folder_location = "./P1"
#the following line will navigate python to the correct foler. chdir stands for change directory.
os.chdir(folder_location)
#the following line return a list of file names in the folder
files = os.listdir()
Now is a good time to notice that we don't actually need the entire image to do this computation. The only thing we really need is the region the particle travels within, as shown below
This is commonly called the Region of Interest (ROI) of an image and we will crop every single image as we import it. Lets do that by importing a test image and checking how we want to crop it:
In [ ]:
test_image = io.imread(files[0])
plt.imshow(test_image['insert slice location'], cmap = cm.gray)
Now that we know where we want to crop the images, we can do so automatically for all of the images in the folder
In [33]:
# the following line creates an empty list that we will populate with the images as we import and crop them.
#That is, ROI is a list of matrices, each one representing an image.
ROIs = []
# the following for-loop imports the images
for image in files:
ROIs.append(io.imread(image)[210::])
# to make sure we are doing things correctly, lets see what the 17th image in the list is
plt.imshow(ROIs[16], cmap=cm.gray)
Out[33]:
Looks good. Okay, now that we have all of the images imported, we want to run a blob finding algorithm that finds the position of the blob in each of the images. However, there are "specles" in the background, so make sure that your blob finding parameters are set correctly.
The goal is for you to use the examples we have provided above to write your own code to systematically find the particles in the ROI list. A general outline of the code is provided below, but if you need further help feel free to ask.
general outline:
In [ ]:
#student solution goes here
We should now have a list of particle locations and sizes. Note that the size information is not important to this application, so we can get rid of it. In general, we find it easier to work with big complicated lists using numpy, so first we will convert the list into a numpy array, and then clean it up.
In [45]:
#the following line converts the list into a numpy array
particles = np.array(particles)
#the following line cleans up array to make it look nicer
particles = particles[::,0]
#this shows us what they array looks like.
print(particles)
Great, using a ROI made this computation faster, and it made sure that the particle was easy to find. Unfortunately it also introduced a small offset error in the verticle position of the particle. This can be seen by taking a look at the picture below:
To correct for this error we simply have to add an offset to the y position of the particle list.
In [47]:
for n, blob in enumerate(particles):
particles[n] = np.add(blob, [210,0,0])
print(particles)
The final step is to turn these pixel locations into measurements of distance between the particle and the center of the rod. To do this, we will use the usual distance formula:
where (y_2, x_2) is the location of the center of the rod and (y_1, y_2) is the location of the particle. From previous measurements, the experimenters know that the center of the rod is approximately at
y_2 = 83
x_2 = 137
so we can change all of the particle location data to particle distance data as follows:
In [80]:
distances = []
for particle in particles:
distance = np.sqrt((83 - particle[0])**2 + (137-particle[1])**2)
distances.append(distance)
distances = np.array(distances)
#print(distances)
Great, so now we have the the distances of the particles from the center of the rod for each movie. The experimenters also know two important pieces of information.
we can use this information to make a plot of the particle's distance as a function of time with proper units.
In [68]:
time = np.linspace(0,len(distances)/2360, len(distances) )
distances = distances*.0059
plt.plot(time,distances)
plt.xlabel('time(s)')
plt.ylabel('distance(mm)')
Out[68]:
You just did particle tracking. Now we'll quickly demonstrate how to find the velocity of the particle.
Recall that the definition of velocity is simply distance/time. Since we know that the time between pictures is 1/2360 seconds, all we have to do is calculate the distance the particle moved between each frame.
In [76]:
velocities = []
for n in range(len(distances)):
if n < (len(distances)-1):
velocity = (distances[n+1] - distances[n])*2360
velocities.append(velocity)
#print(velocities)
Sometimes is it useful/interesting to see the velocity data in a "phase diagram", which is just a plot of the position vs the velocity:
In [79]:
plt.scatter(distances[:-1:],velocities, )
plt.xlabel('distance(mm)')
plt.ylabel('velocity(mm/s)')
Out[79]:
As you can see, the stable particle makes a circle in these "phase diagrams". As you can try (below), the unstable particle will produce sometehing that looks like scribles instead. Phase diagrams are often used to quickly check the stability of a particle, without having to watch the full movie.
In [84]:
os.chdir('..')
folder_location = "./PC"
#the following line will navigate python to the correct foler. chdir stands for change directory.
os.chdir(folder_location)
#the following line return a list of file names in the folder
files_unstable = os.listdir()
#print(files_unstable)