📚 The CoCalc Library - books, templates and other resources
License: OTHER
Introduction to three-dimensional image processing
Images are represented as numpy arrays. A single-channel, or grayscale, image is a 2D matrix of pixel intensities of shape (row, column). We can construct a 3D volume as a series of 2D planes, giving 3D images the shape (plane, row, column). Multichannel data adds a channel dimension in the final position containing color information.
These conventions are summarized in the table below:
| Image type | Coordinates |
|---|---|
| 2D grayscale | (row, column) |
| 2D multichannel | (row, column, channel) |
| 3D grayscale | (plane, row, column) |
| 3D multichannel | (plane, row, column, channel) |
Some 3D images are constructed with equal resolution in each dimension; e.g., a computer generated rendering of a sphere. Most experimental data captures one dimension at a lower resolution than the other two; e.g., photographing thin slices to approximate a 3D structure as a stack of 2D images. The distance between pixels in each dimension, called spacing, is encoded in a tuple and is accepted as a parameter by some skimage functions and can be used to adjust contributions to filters.
Input/Output and display
Three dimensional data can be loaded with skimage.io.imread. The data for this tutorial was provided by the Allen Institute for Cell Science. It has been downsampled by a factor of 4 in the row and column dimensions to reduce computational time.
The distance between pixels was reported by the microscope used to image the cells. This spacing information will be used to adjust contributions to filters and helps decide when to apply operations planewise. We've chosen to normalize it to 1.0 in the row and column dimensions.
To illustrate (no need to read the following cell; execute to generate illustration).
Back to our original data, let's try visualizing the image with skimage.io.imshow.
skimage.io.imshow can only display grayscale and RGB(A) 2D images. We can use skimage.io.imshow to visualize 2D planes. By fixing one axis, we can observe three different views of the image.
Three-dimensional images can be viewed as a series of two-dimensional functions. The display helper function displays 30 planes of the provided image. By default, every other plane is displayed.
Exposure
skimage.exposure contains a number of functions for adjusting image contrast. These functions operate on pixel values. Generally, image dimensionality or pixel spacing does not need to be considered.
Gamma correction, also known as Power Law Transform, brightens or darkens an image. The function is applied to each pixel in the image. A gamma < 1 will brighten an image, while a gamma > 1 will darken an image.
Histogram equalization improves contrast in an image by redistributing pixel intensities. The most common pixel intensities are spread out, allowing areas of lower local contrast to gain a higher contrast. This may enhance background noise.
Most experimental images are affected by salt and pepper noise. A few bright artifacts can decrease the relative intensity of the pixels of interest. A simple way to improve contrast is to clip the pixel values on the lowest and highest extremes. Clipping the darkest and brightest 0.5% of pixels will increase the overall contrast of the image.
Edge detection
Edge detection highlights regions in the image where a sharp change in contrast occurs. The intensity of an edge corresponds to the steepness of the transition from one intensity to another. A gradual shift from bright to dark intensity results in a dim edge. An abrupt shift results in a bright edge.
The Sobel operator is an edge detection algorithm which approximates the gradient of the image intensity, and is fast to compute. skimage.filters.sobel has not been adapted for 3D images. It can be applied planewise to approximate a 3D result.
Filters
In addition to edge detection, skimage.filters provides functions for filtering and thresholding images.
Gaussian filter applies a Gaussian function to an image, creating a smoothing effect. skimage.filters.gaussian takes as input sigma which can be a scalar or a sequence of scalar. This sigma determines the standard deviation of the Gaussian along each axis. When the resolution in the plane dimension is much worse than the row and column dimensions, dividing base_sigma by the image spacing will balance the contribution to the filter along each axis.
Median filter is a noise removal filter. It is particularly effective against salt and pepper noise. An additional feature of the median filter is its ability to preserve edges. This is helpful in segmentation because the original shape of regions of interest will be preserved.
skimage.filters.median does not support three-dimensional images and needs to be applied planewise.
A bilateral filter is another edge-preserving, denoising filter. Each pixel is assigned a weighted average based on neighboring pixels. The weight is determined by spatial and radiometric similarity (e.g., distance between two colors).
skimage.restoration.denoise_bilateral requires a multichannel parameter. This determines whether the last axis of the image is to be interpreted as multiple channels or another spatial dimension. While the function does not yet support 3D data, the multichannel parameter will help distinguish multichannel 2D data from grayscale 3D data.
Thresholding is used to create binary images. A threshold value determines the intensity value separating foreground pixels from background pixels. Foregound pixels are pixels brighter than the threshold value, background pixels are darker. Thresholding is a form of image segmentation.
Different thresholding algorithms produce different results. Otsu's method and Li's minimum cross entropy threshold are two common algorithms. The example below demonstrates how a small difference in the threshold value can visibly alter the binarized image.
Morphological operations
Mathematical morphology operations and structuring elements are defined in skimage.morphology. Structuring elements are shapes which define areas over which an operation is applied. The response to the filter indicates how well the neighborhood corresponds to the structuring element's shape.
There are a number of two and three dimensional structuring elements defined in skimage.morphology. Not all 2D structuring element have a 3D counterpart. The simplest and most commonly used structuring elements are the disk/ball and square/cube.
The most basic mathematical morphology operations are dilation and erosion. Dilation enlarges bright regions and shrinks dark regions. Erosion shrinks bright regions and enlarges dark regions. Other morphological operations are composed of dilation and erosion.
The closing of an image is defined as a dilation followed by an erosion. Closing can remove small dark spots (i.e. “pepper”) and connect small bright cracks. This tends to “close” up (dark) gaps between (bright) features. Morphological opening on an image is defined as an erosion followed by a dilation. Opening can remove small bright spots (i.e. “salt”) and connect small dark cracks. This tends to “open” up (dark) gaps between (bright) features.
These operations in skimage.morphology are compatible with 3D images and structuring elements. A 2D structuring element cannot be applied to a 3D image, nor can a 3D structuring element be applied to a 2D image.
These four operations (closing, dilation, erosion, opening) have binary counterparts which are faster to compute than the grayscale algorithms.
Morphology operations can be chained together to denoise an image. For example, a closing applied to an opening can remove salt and pepper noise from an image.
Functions operating on connected components can remove small undesired elements while preserving larger shapes.
skimage.morphology.remove_small_holes fills holes and skimage.morphology.remove_small_objects removes bright regions. Both functions accept a min_size parameter, which is the minimum size (in pixels) of accepted holes or objects. The min_size can be approximated by a cube.
Segmentation
Image segmentation partitions images into regions of interest. Interger labels are assigned to each region to distinguish regions of interest.
Connected components of the binary image are assigned the same label via skimage.measure.label. Tightly packed cells connected in the binary image are assigned the same label.
A better segmentation would assign different labels to disjoint regions in the original image.
Watershed segmentation can distinguish touching objects. Markers are placed at local minima and expanded outward until there is a collision with markers from another region. The inverse intensity image transforms bright cell regions into basins which should be filled.
In declumping, markers are generated from the distance function. Points furthest from an edge have the highest intensity and should be identified as markers using skimage.feature.peak_local_max. Regions with pinch points should be assigned multiple markers.
After watershed, we have better disambiguation between internal cells.
When cells simultaneous touch the border of the image, they may be assigned the same label. In pre-processing, we typically remove these cells.
Note: This is 3D data---you may not always be able to see connections in 2D!
The watershed algorithm falsely detected subregions in a few cells. This is referred to as oversegmentation.
Plotting the markers on the distance image reveals the reason for oversegmentation. Cells with multiple markers will be assigned multiple labels, and oversegmented. It can be observed that cells with a uniformly increasing distance map are assigned a single marker near their center. Cells with uneven distance maps are assigned multiple markers, indicating the presence of multiple local maxima.
Feature extraction
Feature extraction reduces data required to describe an image or objects by measuring informative features. These include features such as area or volume, bounding boxes, and intensity statistics.
Before measuring objects, it helps to clear objects from the image border. Measurements should only be collected for objects entirely contained in the image.
After clearing the border, the object labels are no longer sequentially increasing. The labels can be renumbered such that there are no jumps in the list of image labels.
skimage.measure.regionprops automatically measures many labeled image features. Optionally, an intensity_image can be supplied and intensity features are extracted per object. It's good practice to make measurements on the original image.
Not all properties are supported for 3D data. Below are lists of supported and unsupported 3D measurements.
skimage.measure.regionprops ignores the 0 label, which represents the background.
Collected measurements can be further reduced by computing per-image statistics such as total, minimum, maximum, mean, and standard deviation.
Perimeter measurements are not computed for 3D objects. The 3D extension of perimeter is surface area. We can measure the surface of an object by generating a surface mesh with skimage.measure.marching_cubes and computing the surface area of the mesh with skimage.measure.mesh_surface_area.
The volume can be visualized using the mesh vertexes and faces.
Challenge problems
Put your 3D image processing skills to the test by working through these challenge problems.
Improve the segmentation
A few objects were oversegmented in the declumping step. Try to improve the segmentation and assign each object a single, unique label. You can try:
generating a smoother image by modifying the
win_sizeparameter inskimage.restoration.denoise_bilateral, or try another filter. Many filters are available inskimage.filtersandskimage.filters.rank.adjusting the threshold value by trying another threshold algorithm such as
skimage.filters.threshold_otsuor entering one manually.generating different markers by changing the size of the footprint passed to
skimage.feature.peak_local_max. Alternatively, try another metric for placing markers or limit the planes on which markers can be placed.
Segment cell membranes
Try segmenting the accompanying membrane channel. In the membrane image, the membrane walls are the bright web-like regions. This channel is difficult due to a high amount of noise in the image. Additionally, it can be hard to determine where the membrane ends in the image (it's not the first and last planes).
Below is a 2D segmentation of the membrane:
The membrane image can be loaded using skimage.io.imread("../images/cells_membrane.tif").
Hint: there should only be one nucleus per membrane.