Image Processing
If you run this python script in a file with *.tiff images, it will
- Convert images to grayscale and invert them (so particles are light instead of dark)
- Subtract a common background computed using the median of the images
- Smooth the images to remove artifacts and accentuate particles
- Apply an intensity threshold to to background subtracted, smoothed image to identify particles
- 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 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)