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)