In [1]:
# import the necessary modules
import numpy as np 
import matplotlib.pyplot as plt
%matplotlib inline

# for pretty plots
import seaborn as sns
sns.set()

# for reading in filenames
import glob

# for importing baterical images
import skimage.io

Let's look a sample image

In [2]:
# read in an example fluor image
im = skimage.io.imread("./data/ecoli_growth/ecoli_TRITC_05.tif")

# take a look at it
im
Out[2]:
array([[205, 205, 205, ..., 205, 205, 206],
       [204, 204, 205, ..., 204, 205, 205],
       [204, 204, 205, ..., 204, 203, 203],
       ...,
       [205, 204, 204, ..., 203, 203, 203],
       [205, 204, 204, ..., 205, 205, 205],
       [205, 204, 204, ..., 205, 205, 205]], dtype=uint16)
In [3]:
# show the image
plt.imshow(im)
plt.grid(False)
In [4]:
# make the image more interactive
%matplotlib notebook

plt.imshow(im)
plt.grid(False)
In [5]:
%matplotlib inline

# generate a histogram of pixel intensities
counts, bins = skimage.exposure.histogram(im)

# plot the histogram
plt.bar(bins, counts)
plt.yscale("log")

# label axes
plt.xlabel("pixel intensity")
plt.ylabel("number of pixels")
Out[5]:
Text(0,0.5,'number of pixels')

Threshold the image

In [6]:
# set the threshold value
threshold = 209

# apply this threshold and show the segmented image
im_thresh = im > threshold 
plt.imshow(im_thresh)
plt.grid(False)

segment all the images

In [7]:
# extract the names of all the phase images
im_names = np.sort(glob.glob('./data/ecoli_growth/ecoli_TRITC_*.tif'))
In [8]:
im_names
Out[8]:
array(['./data/ecoli_growth/ecoli_TRITC_00.tif',
       './data/ecoli_growth/ecoli_TRITC_01.tif',
       './data/ecoli_growth/ecoli_TRITC_02.tif',
       './data/ecoli_growth/ecoli_TRITC_03.tif',
       './data/ecoli_growth/ecoli_TRITC_04.tif',
       './data/ecoli_growth/ecoli_TRITC_05.tif',
       './data/ecoli_growth/ecoli_TRITC_06.tif',
       './data/ecoli_growth/ecoli_TRITC_07.tif',
       './data/ecoli_growth/ecoli_TRITC_08.tif',
       './data/ecoli_growth/ecoli_TRITC_09.tif',
       './data/ecoli_growth/ecoli_TRITC_10.tif',
       './data/ecoli_growth/ecoli_TRITC_11.tif',
       './data/ecoli_growth/ecoli_TRITC_12.tif',
       './data/ecoli_growth/ecoli_TRITC_13.tif',
       './data/ecoli_growth/ecoli_TRITC_14.tif',
       './data/ecoli_growth/ecoli_TRITC_15.tif',
       './data/ecoli_growth/ecoli_TRITC_16.tif',
       './data/ecoli_growth/ecoli_TRITC_17.tif',
       './data/ecoli_growth/ecoli_TRITC_18.tif',
       './data/ecoli_growth/ecoli_TRITC_19.tif',
       './data/ecoli_growth/ecoli_TRITC_20.tif'], dtype='<U38')
In [9]:
# determine the number of frames
n_frames = len(im_names)

# initialize an array for storing bacterial areas
area_tot = np.zeros(n_frames)

# loop through all the images
for n in range(n_frames): 
    
    # read in the frame
    im = skimage.io.imread(im_names[n])
    
    # threshold the image
    im_thresh = im > 209
    
    # save the area in our array
    area_tot[n] = im_thresh.sum()
In [10]:
area_tot
Out[10]:
array([  9269.,  10697.,  12331.,  14680.,  17128.,  19374.,  22730.,
        26903.,  31191.,  34171.,  40581.,  48675.,  58842.,  66354.,
        72622.,  90174., 105100., 127476., 136723., 175299., 197324.])

Plot the area over time and compare to various growth rates

In [11]:
# specify time resolution of imaging
dt = 5 # mins

# get x-axis for time values
times = dt * np.arange(n_frames)

# plot
plt.plot(times,area_tot, '.')

# plot exponential curves on top 
N_0 = area_tot[0]
plt.plot(times,N_0*np.exp(0.025*times))
plt.plot(times,N_0*np.exp(0.03*times))
plt.plot(times,N_0*np.exp(0.035*times))
plt.legend(["data", "k=0.025", "k=0.03", "k=0.035"])
Out[11]:
<matplotlib.legend.Legend at 0x1c21e2c1d0>

try different values of k and see what the corresponding errors are

In [20]:
# the number of differnet k's we want to try
n_points = 1000

# make the range of k's to test
ks = np.linspace(0.01,0.05, n_points)

# array of errors for the differnt values of k
errors = np.zeros(n_points)

# loop through each value of k, compute error
for i in range(n_points):
    # compute theory curve for given k
    theory_data = N_0*np.exp(ks[i]*times)
    
    # compute the error
    error = np.sum((theory_data - area_tot)**2)
    
    # update error array
    errors[i] = error

plot these errors over values of k

In [21]:
plt.plot(ks,errors)
plt.yscale("log")

find the best value of k

In [22]:
k_optimal = np.where(errors == np.min(errors))

k_fit = ks[k_optimal]

print(k_fit)

[0.03054054]
In [23]:
k_optimal
Out[23]:
(array([513]),)