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 to background subtracted, smoothed image to identify particles
5. Write out particle locations (in pixel units) and particle size (# of pixels) for each particle detected

There is an initial calibration factor "cal", which you'll want to set to your calibration in order to have plots calibrated in microns. Set "cal=1" for plotting in pixel units. You may find parts of this code useful for making your own plots.

Another adjustable parameter is "threshold", which determines the intensity level above which a pixel is considered part of a particle. Consult the output pixel_histogram.png to check that this is sensibly set.

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)
#Threshold for intensity when counting a pixel as part of a particle
#From pixel_histogram.png it looks like real particles have intensity >40
threshold=40 

#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
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)