Image Processing

If you run this python script in a file with *.tiff images, it will

1. Convert images to grayscale and invert them (so particles are light instead of dark)
2. Subtract a common background computed using the median of the images
3. Smooth the images to remove artifacts and accentuate particles
4. Apply an intensity threshold to identify particles
5. Write out particle locations (in pixel units) and particle size (# of pixels) for each particle detected

You may find parts of this code useful for doing your own plots.

import numpy as np
import pylab as plt
from scipy import misc, ndimage
from glob import glob

cal=0.167 # calibration microns per pixel (used for images, not *.txt output)

#Read in image
#example = '20140130_224916'
first_tiff = glob('*.tiff')[0] # Get first file
example = first_tiff.split('.')[0] # Remove tiff
bmi = misc.imread('%s.tiff' % example)
#print type(bmi)
print bmi.shape, bmi.dtype

#Convert to grayscale
def rgb2gray(rgb): #This function converts rgb images to grayscale
    r, g, b = rgb[:,:,0], rgb[:,:,1], rgb[:,:,2]
    gray = 0.2989 * r + 0.5870 * g + 0.1140 * b
    return gray
    
bmi = -rgb2gray(bmi) #Negative inverts image
#print bmi.shape 
ny = bmi.shape[0]
nx = bmi.shape[1]

#Plot and save image
plt.imshow(bmi,cmap=plt.cm.gray, extent=[0,nx*cal,0,ny*cal], aspect='auto' )
plt.ylabel(r'y ($\mu$m)')
plt.xlabel(r'x ($\mu$m)')
plt.savefig('%s_grayscale_inverted.png' % example)
plt.clf() #clear figure

#read in all the images
filelist = glob('*.tiff') #Finds all files with .tiff in working directory
#print filelist
bmis=[] #list of images
for filename in filelist:
    im=misc.imread(filename)
    bmis.append(-rgb2gray(im))
    plt.imshow(bmis[-1],cmap=plt.cm.gray, extent=[0,nx*cal,0,ny*cal], aspect='auto' )
    plt.ylabel(r'y ($\mu$m)')
    plt.xlabel(r'x ($\mu$m)')
    plt.savefig('%s_grayscale_inverted.png' % (filename.split('.')[0]))


#Compute background by finding the median intensity value for each pixel
background = np.median(bmis,axis=0)
plt.imshow(background,cmap=plt.cm.gray, extent=[0,nx*cal,0,ny*cal], aspect='auto' )
plt.ylabel(r'y ($\mu$m)')
plt.xlabel(r'x ($\mu$m)')
plt.savefig('background.png')
plt.clf()

#Subtract background
bmi=bmi-background
plt.imshow(bmi,cmap=plt.cm.gray, extent=[0,nx*cal,0,ny*cal], aspect='auto' )
plt.ylabel(r'y ($\mu$m)')
plt.xlabel(r'x ($\mu$m)')
plt.savefig('%s_grayscale_inverted_bgsubtracted.png' % example)
plt.clf()
for i in range(len(bmis)):
    bmis[i] = bmis[i]-background

#Smooth to particle size
bmi = ndimage.uniform_filter(bmi, size=7)
plt.imshow(bmi,cmap=plt.cm.gray, extent=[0,nx*cal,0,ny*cal], aspect='auto' )
plt.ylabel(r'y ($\mu$m)')
plt.xlabel(r'x ($\mu$m)')
plt.savefig('%s_grayscale_inverted_bgsubtracted_smoothed.png' % example)
plt.clf()
for i in range(len(bmis)):
    bmis[i] = ndimage.uniform_filter(bmis[i], size=5)


#Choose threshold for particles
plt.hist(bmi.ravel(), bins=100)
plt.xlabel('Intensity')
plt.ylabel('# of Pixels')
ax=plt.gca()
ax.set_yscale('log')
plt.savefig('pixel_histogram.png')
plt.clf()

#From histogram looks like real particles have brightness higher than 40
#Create binary images: '1' if intensity is greater than 40, '0' otherwise
threshold=40
inds = np.where(bmi>threshold)
bmi = np.zeros(bmi.shape)
bmi[inds]=1.
plt.imshow(bmi,cmap=plt.cm.gray, extent=[0,nx*cal,0,ny*cal], aspect='auto' )
plt.ylabel(r'y ($\mu$m)')
plt.xlabel(r'x ($\mu$m)')
plt.savefig('%s_grayscale_inverted_bgsubtracted_smoothed_threshold.png' % example)
plt.clf()
for i in range(len(bmis)):
    inds = np.where(bmis[i]>threshold)
    bmis[i] = np.zeros(bmi.shape)
    bmis[i][inds] = 1.
    plt.imshow(bmis[i],cmap=plt.cm.gray, extent=[0,nx*cal,0,ny*cal], aspect='auto' )
    plt.ylabel(r'y ($\mu$m)')
    plt.xlabel(r'x ($\mu$m)')
    plt.savefig('%s_grayscale_inverted_bgsubtracted_smoothed_threshold.png' % (filelist[i].split('.')[0]))
    plt.clf()

#Find particles

#The label function find groups of contiguous pixels with '1' values
#It returns an array 'labels', in which the '1' values are replaced with
#distinct integers for each group
#
# Example, lets say there is a 1D array: '0,1,1,0,1,0,1,1,1'
# The labels output would be '0,1,1,0,2,0,3,3,3'
#

labels, num = ndimage.measurements.label(bmi)
print "Found %d particles." % num
for i in range(num):
    inds = np.where(labels==i+1) #get indices corresponding to (i+1)th particle
    print i, inds[0].mean(), inds[1].mean(), len(inds[0])

for i in range(len(bmis)):
    labels, num = ndimage.measurements.label(bmis[i])
    print "Found %d particles in %s" % (num, filelist[i])
    particles = []
    for j in range(num):
        inds = np.where(labels==j+1)
        particles.append([inds[0].mean(), inds[1].mean(), len(inds[0])])
    np.savetxt("%s.txt" % (filelist[i].split('.')[0]), particles)