Exercise: a growth curve from real images

Exercise: a growth curve from real images#

Duration ~40 min · Session Day 2, Python notebooks

This is the whole course in one exercise: segment a series of images, count the objects, fit a model to the counts, and report a biological quantity.

The data is real, and it does not behave as neatly as the simulated data in the walkthrough.

Data: data/misc/timelapse/: seven phase-contrast frames, 30 minutes apart.

import numpy as np
import tifffile
import matplotlib.pyplot as plt

from skimage.measure import label
from skimage.color import label2rgb
from skimage.morphology import remove_small_objects, closing, disk
from scipy.ndimage import binary_fill_holes
from scipy.optimize import curve_fit
import os; os.environ["DISPLAY"] = ":1"

from course import DATA
from iaf.plot import imshow, show_labels

frames = sorted((DATA / "misc" / "timelapse").glob("*.tif"))
print(len(frames), "frames")

fig, axes = plt.subplots(1, 7, figsize=(21, 3.2))
for ax, path in zip(axes, frames):
    imshow(tifffile.imread(path), title=path.stem[-4:], ax=ax, auto_stretch=True)
plt.show()

Task 1: why the usual threshold will not work#

Print the min, max, mean and standard deviation of the first and last frames, then try an Otsu threshold on each and count the objects.

Hint

threshold_otsu from skimage.filters, then label(image > t).max().

from skimage.filters import threshold_otsu

# --- your turn ---
for index in (0, 6):
    image = tifffile.imread(frames[index]).astype(float)
    t = ...        # TODO: Otsu threshold
    mask = ...     # TODO
    print(...)     # TODO: report min, max, mean, sd, % foreground, object count

Otsu always splits the histogram into two groups whatever it is given. Look at the standard deviations. Why does it fail here, and what is different about phase-contrast images?

(edit this cell to add your answer)

If you are stuck

In fluorescence, objects are bright and background is dark, so one threshold separates them. In phase contrast the background is a uniform mid-gray and each cell has a dark body with a bright halo, so cells are both brighter and darker than the background, and no single cut-off can capture them.

Task 2: segment on deviation instead#

Since cells deviate from the background in both directions, threshold the absolute difference from the background level rather than the raw intensity.

Write a count_cells(image) function and check it on the first and last frames.

Hint 1: the background level

The background dominates the image, so np.median(image) is a robust estimate of it.

Hint 2: the binary image

np.abs(image - background) > 4. Then tidy it up: closing(mask, disk(2)) to join a cell body to its halo, binary_fill_holes, and remove_small_objects.

# --- your turn ---
def count_cells(image, deviation=4, min_size=30):
    background = ...   # TODO: robust background level
    mask = ...         # TODO: pixels deviating from it in either direction
    mask = ...         # TODO: closing, fill holes, remove small objects
    return label(mask).max()

for index in (0, 6):
    image = tifffile.imread(frames[index]).astype(float)
    print(f"frame {index + 1}: {count_cells(image)} cells")

Task 3: check the parameter is not doing the work#

Before trusting the counts, vary the deviation cut-off and confirm the answer does not depend on it. A result that swings with an arbitrary parameter is not a measurement. Visualize the masks you found to identify the appropriate deviation parameters

Hint

Loop over deviation in (0.5, 1, 1.5, 2, 2.5, 3, 3.5) and print the count for each frame. Additionally overlay the mask with the raw image of timepoint 4

# --- your turn ---
deviations = (0.5, 1, 1.5, 2, 2.5, 3, 3.5)
lbl = []
img3_counts = []
images = [tifffile.imread(path).astype(float) for path in frames]

for deviation in deviations:
    counts = ...   # TODO: count every frame at this cut-off
    print(...)

At this point you should have a table summarizing the counts for all image frames and thresholds.

How much do the counts move as the cut-off changes?
Would you be comfortable reporting these numbers in a publication?

(edit this cell to add your answer)

Task 4: fit the growth curve#

Count every frame, and fit an exponential. From the microscopist who took the images, you know the image frames are 30 minutes apart.

Hint 1: the model

def exponential(t, n0, k):
    return n0 * np.exp(k * t)

Hint 2: the fit

curve_fit(exponential, time, counts, p0=(1, 0.02)), then take the parameter uncertainties from np.sqrt(np.diag(covariance)).

# --- your turn ---
counts = ...   # TODO: count every frame
time = ...     # TODO: 30 minutes per frame

def exponential(t, n0, k):
    return ...

parameters, covariance = ...   # TODO
print(f"doubling time = {...:.1f} min")

Task 5: check the residuals#

Plot the data with the fitted curve, and the residuals beside it.

Hint

Two panels, as in section 6 of the walkthrough. Residual = counts - exponential(time, *parameters).

# --- your turn ---
... two panels: data with fit, and residuals ...

Try to interpret the plot.

Are the residuals random, or is there structure?
What does that tell you about the model?

(edit this cell to add your answer)

Answer:

Structure, not random: The residuals show a distinct pattern (four negative values, followed by two strong positive values, ending with a negative drop). They are not randomly distributed around zero.
What it tells you: The model is misspecified and systematically under- or overestimating specific segments of your data. This indicates that your current model is failing to capture the true underlying trend, likely requiring a non-linear approach (such as adding a quadratic term) or a missing variable.

Task 6: the same model, fitted two ways#

There is a second way to fit an exponential. Taking logs turns it into a straight line,

\[\ln N = \ln N_0 + kt\]

so you can fit it with np.polyfit instead. Do that, and compare the doubling time with what you got in task 4.

Hint

slope, intercept = np.polyfit(time, np.log(counts), 1): note polyfit returns the highest-order coefficient first.

# --- your turn ---
slope, intercept = ...   # TODO: fit log(counts) against time

print(f"fit in linear space : doubling = {...:.1f} min")
print(f"fit in log space    : doubling = {...:.1f} min")

You found that the two doubling times differ by several minutes, from the same data and the same model.

Which would you report?

Look at the two panels before answering. On the linear scale, squared error is dominated by the last few points, because they are the largest, so that fit prioritises the end of the curve. On the log scale every point contributes equally in relative terms, so the early points matter as much as the late ones.

(edit this cell to add your answer)

Answer: Neither is wrong. They are answers to slightly different questions, and the disagreement is only large because the exponential model does not fit this data well in the first place.

Task 7: a model that fits#

Growth cannot continue forever; a dish runs out of room. The logistic model adds a carrying capacity \(K\):

\[N(t) = \frac{K}{1 + \frac{K - N_0}{N_0} e^{-kt}}\]

Fit it, and compare the residuals with the exponential’s.

Hint

Supply a sensible p0: the carrying capacity is at least as large as your final count, n0 is around 1, and k is roughly what you already fitted.

# --- your turn ---
def logistic(t, carrying_capacity, n0, k):
    return ...

logistic_parameters, _ = ...   # TODO, with a sensible p0

... compare R2 and residuals for both models ...

The logistic model fits better.

Does that make it the right choice here? Consider that it has three parameters against the exponential’s two, and that the fitted carrying capacity is an extrapolation well beyond the last frame you actually observed.

(edit this cell to add your answer)

Answer: A more flexible model will almost always fit better. That is not, by itself, a reason to prefer it. The questions worth asking are whether the extra parameter is biologically motivated, a carrying capacity is, and whether it is constrained by data you actually have. Here it is not: the cells are still growing in the last frame, so K is being inferred from curvature alone. Report it with that caveat, or not at all.

To evaluate this choice objectively, we should look at a quantitative metric like the Bayesian Information Criterion (BIC). BIC mathematically penalizes models for adding extra parameters (3 vs. 2) to prevent overfitting. If the logistic model yields a significantly lower BIC despite this penalty, it suggests the improvement in fit justifies the extra complexity. However, if the BIC values are close, the simpler exponential model remains the safer, more parsimonious choice due to the dangerous extrapolation of

Exercise: a growth curve from real images

Contents

Exercise: a growth curve from real images#

Task 1: why the usual threshold will not work#

Task 2: segment on deviation instead#

Task 3: check the parameter is not doing the work#

Task 4: fit the growth curve#

Task 5: check the residuals#

Task 6: the same model, fitted two ways#

Task 7: a model that fits#