>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> import scipy as sp

1.5.11.4. Example of solution for the image processing exercise: unmolten grains in glass

1. Open the image file MV_HFV_012.jpg and display it. Browse through the keyword arguments in the docstring of imshow to display the image with the “right” orientation (origin in the bottom left corner, and not the upper left corner as for standard arrays).

   >>> dat = plt.imread('data/MV_HFV_012.jpg')
   

2. Crop the image to remove the lower panel with measure information.

   >>> dat = dat[:-60]
   

3. Slightly filter the image with a median filter in order to refine its histogram. Check how the histogram changes.

   >>> filtdat = sp.ndimage.median_filter(dat, size=(7,7))
   >>> hi_dat = np.histogram(dat, bins=np.arange(256))
   >>> hi_filtdat = np.histogram(filtdat, bins=np.arange(256))
   

   

4. Using the histogram of the filtered image, determine thresholds that allow to define masks for sand pixels, glass pixels and bubble pixels. Other option (homework): write a function that determines automatically the thresholds from the minima of the histogram.

   >>> void = filtdat <= 50
   >>> sand = np.logical_and(filtdat > 50, filtdat <= 114)
   >>> glass = filtdat > 114
   

5. Display an image in which the three phases are colored with three different colors.

   >>> phases = void.astype(int) + 2*glass.astype(int) + 3*sand.astype(int)
   

   

6. Use mathematical morphology to clean the different phases.

   >>> sand_op = sp.ndimage.binary_opening(sand, iterations=2)
   

7. Attribute labels to all bubbles and sand grains, and remove from the sand mask grains that are smaller than 10 pixels. To do so, use sp.ndimage.sum or np.bincount to compute the grain sizes.

   >>> sand_labels, sand_nb = sp.ndimage.label(sand_op)
   >>> sand_areas = np.array(sp.ndimage.sum(sand_op, sand_labels, np.arange(sand_labels.max()+1)))
   >>> mask = sand_areas > 100
   >>> remove_small_sand = mask[sand_labels.ravel()].reshape(sand_labels.shape)
   

   

8. Compute the mean size of bubbles.

   >>> bubbles_labels, bubbles_nb = sp.ndimage.label(void)
   >>> bubbles_areas = np.bincount(bubbles_labels.ravel())[1:]
   >>> mean_bubble_size = bubbles_areas.mean()
   >>> median_bubble_size = np.median(bubbles_areas)
   >>> mean_bubble_size, median_bubble_size
   (np.float64(1699.875), np.float64(65.0))