The aim of this IPython Notebook is to introduce you to how real quantitative analysis of an image is performed. In order to do this, we will complete the following analysis procedure:
We are going to start by loading a raw image file that is representative of a data set that may have been acquired during an imaging experiment. To do this, we will use the Scipy (http://www.scipy.org/) library.
In general, image data can be stored in either raster or vector graphics formats. Vector graphics scale with the resolution of the device on which they are viewed, whereas the appearance of raster graphics is resolution dependent and cannot scale arbitrarily. Vector graphics are typically used for graphical design applications whilst raster graphics are used for photographic images.
In this demonstration we will be concerned with raster graphics which are effectively a dot matrix of data. As medical images may be optical or non-optical (e.g. X-ray) we will also discuss the coversion between colour and grayscale. At each point in the matrix of a colour (RGB) raster image a red (R), green (G) and blue (B) value is stored.
The raster image we will load is of Sally Ride (http://en.wikipedia.org/wiki/Sally_Ride) the first US female in space. We will use this image as our phantom (test data).
In [1]:
# Import required libraries.
# Get the Scipy ndimage library.
from scipy import ndimage
# Specify the path to the image (PNG) file.
sallyPath = "./Sally.png"
# Load the image file into memory.
sallyPNG = ndimage.imread(sallyPath)
Now we have loaded the two data sets, we can investigate their properties.
In [2]:
# Check the shape of the image object
print "The shape of the external file (PNG) structure is %d x %d x %d." %(sallyPNG.shape[0], sallyPNG.shape[1], sallyPNG.shape[2])
print "The maximum pixel intensity is %d." %(sallyPNG.max())
print "The data type is %s." %(sallyPNG.dtype)
We can see that the PNG file has three values at each location in the 512 x 512 matrix, corresponding to the intensity of R, G and B. As a side note, an RGB image may have three or four columns. The fourth column would correspond to an $\alpha$-value (A) which describes the opacity of the image and is utilised for overlaying multiple images. We will not worry about $\alpha$-values in this demonstration. A PNG would typically be 32-bit (four bytes) per pixel for RGBA, with 8-bit integer precision per channel (R, G, B and A), per pixel. For grayscale, the PNG might have 8-bit integer precision per pixel with a byte-depth of only one.
Now we will plot the image using the imshow function in Matplotlib. As we are plotting data, we must also be concerned with presentation. Matplotlib does not have a layout-engine and so certain formatting tasks must be undertaken by the user. In the following example, the first section of code is present to overwrite the default rcParams dictionary with those suitable for the IPython Notebook.
In [3]:
# Specify that we want to show the plots in the notebook.
%matplotlib inline
# Get the Matplotlib plotting libraries.
import pylab as plt
# Plot the PNG in the jet colourmap.
plt.figure()
plt.imshow(sallyPNG, interpolation = 'nearest', origin='upper', cmap="jet")
plt.xlabel("Pixel number")
plt.ylabel("Pixel number")
plt.clim([0,255])
plt.colorbar()
plt.show()
In the plot above all three colour components (R, G and B) are shown as is typical for a standard display. We have specified the "jet" colour map explicitly in this example, although "jet" is the defaults colour map in Matplotlib.
The image we have imported has three color channels; we can combine these into a single grayscale channel so that we have a monochrome image to work with. The data that we would acquire from most typical experiments is expressed in terms of intensities (grayscale), not RGB triplets.
As a side note, it is useful to understand that grayscale images can be stored as integer (e.g. 8-bit PNG) or floating point (e.g. 32-bit TIFF - complicated container not discussed here). Scipy and Matplotlib both call the Python Image Library (PIL) for image reading and writing. The PIL can only read up to 16-bit TIFFs.
In [4]:
# A more elegant way to separate the image into its RGB components is to use the transpose of the array of input data.
redSallyT, greenSallyT, blueSallyT = ndimage.imread(sallyPath).T
# Now apply step to of the conversion from sRGB to linear grayscale.
graySally = (0.299 * redSallyT) + (0.587 * greenSallyT) + (0.114 * blueSallyT)
# We have to apply a second transpose as the initial transpose swapped the X and Y axes as well as the Z axis.
graySally = graySally.T.astype("uint8")
In [5]:
# Plot the PNG in the jet colourmap.
# Can't keep the same figure size for multi-sub-plot so need to specify something sensible.
plt.figure(figsize = (12.0, 4.0))
plt.subplot(1,2,1)
plt.imshow(sallyPNG, interpolation = 'nearest', origin='upper', cmap="jet")
plt.xlabel("Pixel number")
plt.ylabel("Pixel number")
plt.clim([0,255])
plt.colorbar()
# Plot the grayscale image in grayscale.
plt.subplot(1,2,2)
plt.imshow(graySally, interpolation = 'nearest', origin='upper', cmap="gray")
plt.xlabel("Pixel number")
plt.ylabel("Pixel number")
plt.clim([0,255])
plt.colorbar()
plt.show()
We have However, we have only populated one of the 3 subplots. What we really need is a function that incorporates the syntax we have just used but that populates the 3 subplots using a for loop.
Write a function called plotNby1Subplot that takes a list of N tuples of images with their corresponding colour-maps and plots them in a 1xN row of subplots. Demonstrate that it correctly plots the tuples:
(sallyPNG, "jet"),
(graySally, "jet"),
(graySally, "gray").Use the same plotNby1Subplot function and the correct colour-maps to display the R, G and B components of sallyPNG.
In [6]:
# SOLUTION TO CHALLENGE 1
def plotNby1Subplot(icts, figSize = (12.0, 4.0)):
# Get handles for the figure and the subplots.
fig = plt.figure(figsize = figSize)
for index, ict in enumerate(icts):
# Plot the image.
plt.subplot(1, len(icts), index + 1)
im2Plot = icts[index][0]
plt.imshow(im2Plot, interpolation = 'nearest', origin = 'upper', cmap = icts[index][1])
plt.xlabel("Pixel number")
plt.ylabel("Pixel number")
plt.clim([0,255])
plt.colorbar()
# Define the tuple data.
imageColbarTupList = [(sallyPNG, "jet"),
(graySally, "jet"),
(graySally, "gray")]
# Plot the images
plotNby1Subplot(imageColbarTupList, figSize = (14.0, 3.0))
plt.show()
In [7]:
# SOLUTION TO CHALLENGE 2
# We can access the R component of the array.
redSally = sallyPNG[:,:,0]
# We can also access the G and B components:
greenSally = sallyPNG[:,:,1]
blueSally = sallyPNG[:,:,2]
imageColbarTupList = [(redSally, "Reds"),
(greenSally, "Greens"),
(blueSally, "Blues")]
# Plot the images
plotNby1Subplot(imageColbarTupList, figSize = (14.0, 3.0))
plt.show()
Up to now all that we have done is load an image file and present the data to explain how it is stored in different image files. RGB data is not usually the most useful format for quantitative analysis in medical imaging. What we need are the grayscale intensities, so that we have a single number at each point in the 512 x 512 matrix.
The rigorous way to convert RGB to grayscale follows a multi-step process (http://en.wikipedia.org/wiki/SRGB):
We are not going to worry about the gamma correction and will just apply Step 2 which is a good enough approximation.
In [8]:
def plotSingleFigImshow(im2Plot):
# Plot grayscale image.
plt.imshow(im2Plot, interpolation = 'nearest', origin='upper', cmap='gray')
plt.xlabel("Pixel numer")
plt.ylabel("Pixel numer")
plt.clim([0,255])
plt.colorbar()
# Plot the image.
plotSingleFigImshow(graySally)
plt.show()
Now we have the grayscale (intensity) distribution of the image we are ready to analyse the data.
In image processing, we often wish to apply a point operation to the data set. In this example we will set a window with a level of 50 and a width of 60 ($\pm$30). All pixels that lie outside of the window (20 < I < 80) will be set to zero. We could do this with nested-loops, however, this is cumbersome and messy. A neater way is to use a boolean mask:
In [9]:
# Generate a mask by simple thresholding.
boolMask = (graySally > 80) | (graySally < 20)
# In the mask, all pixels with intensity >200 OR <100 will be set to True.
# Now set these pixels to zero.
graySally[boolMask] = 0
# Plot the image.
plotSingleFigImshow(graySally)
plt.show()
We can see that all regions of the data that lie outside the window now appear black in the image.
If we now decide that we want to apply a different threshold we will run into a problem. We will demonstrate this by applying a window with a level of 120 and a width of 40 ($\pm$20).
In [10]:
# Generate a mask by simple thresholding.
boolMask = (graySally > 140) | (graySally < 100)
# In the mask, all pixels with intensity >200 OR <100 will be set to True.
# Now set these pixels to zero.
graySally[boolMask] = 0
# Plot the image.
plotSingleFigImshow(graySally)
plt.show()
We can see that the image is black. That is because we are applying a second threshold to the same instance of the data that we applied the initial window to. Whenever we want to manipulate data we should first make a copy so that we maintain the integrity of the initial instance. We can do this by using the copy function in the copy library:
In [11]:
# Import the copy library.
import copy as cp
# We need to go back and regenerate graySally.
graySally = (0.299 * redSally) + (0.587 * greenSally) + (0.114 * blueSally)
# And get it back in the right data type
graySally = graySally.astype("uint8")
# We must make a copy of the original image otherwise we will not be able to compare to the original.
copy1Sally = graySally.copy()
copy2Sally = graySally.copy()
# Generate a mask by simple thresholding.
boolMask = (graySally > 80) | (graySally < 20)
# In the mask, all pixels with intensity >80 OR <20 will be set to True.
# Now set these pixels to zero.
copy1Sally[boolMask] = 0
# Generate a mask by simple thresholding.
boolMask = (graySally > 140) | (graySally < 100)
# In the mask, all pixels with intensity >140 OR <100 will be set to True.
# Now set these pixels to zero.
copy2Sally[boolMask] = 0
# Plot the images
# If you managed to write the function in the challenge, use this to plot the figures. Otherwise, revert to the non-ideal method.
try:
# Define the tuple data.
imageColbarTupList = [(copy1Sally, "gray"),
(copy2Sally, "gray")]
plotNby1Subplot(imageColbarTupList)
except:
print "Plotting function does not exist. Using non-ideal method"
plt.figure(figsize = (4.0, 2.0))
# Plot the grayscale image copy1Sally.
plt.subplot(1,3,1)
plt.imshow(sallyPNG, interpolation = 'nearest', origin='upper', cmap="jet")
plt.xlabel("Pixel number")
plt.ylabel("Pixel number")
plt.clim([0,255])
plt.colorbar()
# Plot the grayscale image copy2Sally.
plt.subplot(1,3,2)
plt.imshow(graySally, interpolation = 'nearest', origin='upper', cmap="jet")
plt.xlabel("Pixel number")
plt.ylabel("Pixel number")
plt.clim([0,255])
plt.colorbar()
plt.show()