Welcome to the quickstart guide for dask-image.
We first install the library scikit-image for easier access to the example image data there.
If you are running this notebook on your own computer and not in the mybinder environment, you'll additionally need to ensure your Python environment contains:
You can refer to the full list of dependencies used for the
dask-examples repository, available in the
binder/environment.yml file here (note that the
nomkl package is not available for Windows users): https://github.com/dask/dask-examples/blob/main/binder/environment.yml
When you import dask-image, be sure to use an underscore instead of a dash between the two words.
import dask_image.imread import dask_image.ndfilters import dask_image.ndmeasure import dask.array as da
We'll also use matplotlib to display image results in this notebook.
import matplotlib.pyplot as plt %matplotlib inline
We'll use some example image data from the scikit-image library in this tutorial. These images are very small, but will allow us to demonstrate the functionality of dask-image.
Let's download and save a public domain image of the astronaut Eileen Collins to a temporary directory. This image was originally downloaded from the NASA Great Images database https://flic.kr/p/r9qvLn, but we'll access it with scikit-image's
import os from skimage import data, io output_filename = os.path.join('temp', 'astronaut.png') io.imsave(output_filename, data.astronaut())
Really large datasets often can't fit all of the data into a single file, so we'll chop this image into four parts and save the image tiles to a second temporary directory. This will give you a better idea of how you can use dask-image on a real dataset.
io.imsave(os.path.join('temp-tiles', 'image-00.png'), data.astronaut()[:256, :256, :]) # top left io.imsave(os.path.join('temp-tiles', 'image-01.png'), data.astronaut()[:256, 256:, :]) # top right io.imsave(os.path.join('temp-tiles', 'image-10.png'), data.astronaut()[256:, :256, :]) # bottom left io.imsave(os.path.join('temp-tiles', 'image-11.png'), data.astronaut()[256:, 256:, :]) # bottom right
Now we have some data saved, let's practise reading in files with dask-image and processing our images.
import os filename = os.path.join('temp', 'astronaut.png') print(filename)
astronaut = dask_image.imread.imread(filename) print(astronaut) plt.imshow(astronaut[0, ...]) # display the first (and only) frame of the image
This has created a dask array with
shape=(1, 512, 512, 3). This means it contains one image frame with 512 rows, 512 columns, and 3 color channels.
Since the image is relatively small, it fits entirely within one dask-image chunk, with
chunksize=(1, 512, 512, 3).
In many cases, you may have multiple images stored on disk, for example:
image_NN.png. These can be read into a dask array as multiple image frames.
Here we have the astronaut image split into four non-overlapping tiles:
image_00.png= top left image (index 0,0)
image_01.png= top right image (index 0,1)
image_10.png= bottom left image (index 1,0)
image_11.png= bottom right image (index 1,1)
This filename pattern can be matched with regex:
import os filename_pattern = os.path.join('temp-tiles', 'image-*.png') tiled_astronaut_images = dask_image.imread.imread(filename_pattern) print(tiled_astronaut_images)
This has created a dask array with
shape=(4, 256, 256, 3). This means it contains four image frames; each with 256 rows, 256 columns, and 3 color channels.
There are four chunks in this particular case. Each image frame here is a separate chunk with
chunksize=(1, 256, 256, 3).
fig, ax = plt.subplots(nrows=2, ncols=2) ax[0,0].imshow(tiled_astronaut_images) ax[0,1].imshow(tiled_astronaut_images) ax[1,0].imshow(tiled_astronaut_images) ax[1,1].imshow(tiled_astronaut_images) plt.show()
Next you'll want to do some image processing, and apply a function to your images.
We'll use a very simple example: converting an RGB image to grayscale. But you can also use this method to apply arbittrary functions to dask images. To convert our image to grayscale, we'll use the equation to calculate luminance (reference pdf)":
Y = 0.2125 R + 0.7154 G + 0.0721 B
We'll write the function for this equation as follows:
def grayscale(rgb): result = ((rgb[..., 0] * 0.2125) + (rgb[..., 1] * 0.7154) + (rgb[..., 2] * 0.0721)) return result
Let's apply this function to the astronaut image we read in as a single file and visualize the computation graph.
(Visualizing the computation graph isn't necessary most of the time but it's helpful to know what dask is doing under the hood, and it can also be very useful for debugging problems.)
single_image_result = grayscale(astronaut) print(single_image_result) single_image_result.visualize()
We also see that there are no longer three color channels in the shape of the result, and that the output image is as expected.
print("Original image dimensions: ", astronaut.shape) print("Processed image dimensions:", single_image_result.shape) fig, (ax0, ax1) = plt.subplots(nrows=1, ncols=2) ax0.imshow(astronaut[0, ...]) # display the first (and only) frame of the image ax1.imshow(single_image_result[0, ...], cmap='gray') # display the first (and only) frame of the image # Subplot headings ax0.set_title('Original image') ax1.set_title('Processed image') # Don't display axes ax0.axis('off') ax1.axis('off') # Display images plt.show(fig)
The syntax is identical to apply a function to multiple images or dask chunks. This is an example of an embarrassingly parallel problem, and we see that dask automatically creates a computation graph for each chunk.
result = grayscale(tiled_astronaut_images) print(result) result.visualize()
Let's take a look at the results.
fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(nrows=2, ncols=2) ax0.imshow(result[0, ...], cmap='gray') ax1.imshow(result[1, ...], cmap='gray') ax2.imshow(result[2, ...], cmap='gray') ax3.imshow(result[3, ...], cmap='gray') # Subplot headings ax0.set_title('First chunk') ax1.set_title('Second chunk') ax2.set_title('Thurd chunk') ax3.set_title('Fourth chunk') # Don't display axes ax0.axis('off') ax1.axis('off') ax2.axis('off') ax3.axis('off') # Display images plt.show(fig)
OK, Things are looking pretty good! But how can we join these image chunks together?
So far, we haven't needed any information from neighboring pixels to do our calculations. But there are lots of functions (like those in dask-image ndfilters) that do need this for accurate results. You could end up with unwanted edge effects if you don't tell dask how your images should be joined.
Dask has several ways to join chunks together: Stack, Concatenate, and Block.
Block is very versatile, so we'll use that in this next example. You simply pass in a list (or list of lists) to tell dask the spatial relationship between image chunks.
data = [[result[0, ...], result[1, ...]], [result[2, ...], result[3, ...]]] combined_image = da.block(data) print(combined_image.shape) plt.imshow(combined_image, cmap='gray')
Most analysis pipelines require some degree of image preprocessing. dask-image has a number of inbuilt filters available via dask-image ndfilters
Commonly a guassian filter may used to smooth the image before segmentation. This causes some loss of sharpness in the image, but can improve segmentation quality for methods that rely on image thresholding.
smoothed_image = dask_image.ndfilters.gaussian_filter(combined_image, sigma=[1, 1])
We see a small amount of blur in the smoothed image.
fig, (ax0, ax1) = plt.subplots(nrows=1, ncols=2) ax0.imshow(smoothed_image, cmap='gray') ax1.imshow(smoothed_image - combined_image, cmap='gray') # Subplot headings ax0.set_title('Smoothed image') ax1.set_title('Difference from original') # Don't display axes ax0.axis('off') ax1.axis('off') # Display images plt.show(fig)
Since the gaussian filter uses information from neighbouring pixels, the computational graph looks more complicated than the ones we looked at earlier. This is no longer embarrassingly parallel. Where possible dask keeps the computations for each of the four image chunks separate, but must combine information from different chunks near the edges.
After the image preprocessing, we segment regions of interest from the data. We'll use a simple arbitrary threshold as the cutoff, at 75% of the maximum intensity of the smoothed image.
threshold_value = 0.75 * da.max(smoothed_image).compute() print(threshold_value)
threshold_image = smoothed_image > threshold_value plt.imshow(threshold_image, cmap='gray')
Next, we label each region of connected pixels above the threshold value. For this we use the
label function from dask-image ndmeasure. This will return both the label image, and the number of labels.
label_image, num_labels = dask_image.ndmeasure.label(threshold_image)
print("Number of labels:", int(num_labels)) plt.imshow(label_image, cmap='viridis')
There are a number of inbuilt functions in dask-image ndmeasure useful for quantitative analysis.
We'll use the
dask_image.ndmeasure.standard_deviation() functions, and apply them to each label region with
index = list(range(int(num_labels))) # Note that we're including the background label=0 here, too. out_dtype = float # The data type we want to use for our results. default = None # The value to return if an element of index does not exist in the label image. mean_values = dask_image.ndmeasure.labeled_comprehension(combined_image, label_image, index, dask_image.ndmeasure.mean, out_dtype, default, pass_positions=False) print(mean_values.compute())
Since we're including label 0 in our index, it's not surprising that the first mean value is so much lower than the others - it's the background region below our cutoff threshold for segmentation.
Let's also calculate the standard deviation of the pixel values in our greyscale image.
stdev_values = dask_image.ndmeasure.labeled_comprehension(combined_image, label_image, index, dask_image.ndmeasure.standard_deviation, out_dtype, default, pass_positions=False)
Finally, let's load our analysis results into a pandas table and then save it as a csv file.
import pandas as pd df = pd.DataFrame() df['label'] = index df['mean'] = mean_values.compute() df['standard_deviation'] = stdev_values.compute() df.head()
df.to_csv('example_analysis_results.csv') print('Saved example_analysis_results.csv')
I hope this guide has helped you to get started with dask-image.
The dask-examples repository has a number of other example notebooks: https://github.com/dask/dask-examples
Scaling up with dask distributed
Saving image data with zarr
In some cases it may be necessary to save large data after image processing, zarr is a python library that you may find useful.
You recall we saved some example data to the directories
temp-tiles/. To delete the contents, run the following command:
!rm -r temp
!rm -r temp-tiles