Apply operations on data

One goal behind the clesperanto project is to standardise as much as possible the library usage in order to make is as easy as possible to learn and to use. In this notebook we will see how to apply operations on a data.

All operations follows the same pattern:

• cle.operation_name(input, output, parameters)

where:

• cle is the librairy handle

• operation_name is the name of the operation you want to apply

• input is the input data

• output is the output data

• parameters are the parameters of the operation

Of course, the method signature can becomme more complex if you want to apply more advanced operations, so do not hesitate to check the documentation of the operation you want to use.

    cle.operation_name?

Like always let’s first load an image …

import pyclesperanto_prototype as cle
import numpy as np
from skimage.io import imread, imshow
cle.select_device('TX')  # TODO: change to your GPU

<NVIDIA GeForce RTX 2080 SUPER on Platform: NVIDIA CUDA (2 refs)>

image = imread("../../data/blobs.tif")
imshow(image)

<matplotlib.image.AxesImage at 0x7f2b2180bd50>

Now, to apply a gaussian blur on the image, we will need to first push the image into the GPU memory, create an output memory space in the GPU memory, and apply the operations. After that we only need to pull the image back to the CPU memory in order to display it.

cle.gaussian_blur?

Signature:
cle.gaussian_blur(
    source: Union[numpy.ndarray, pyclesperanto_prototype._tier0._pycl.OCLArray, pyopencl._cl.Image, pyclesperanto_prototype._tier0._pycl._OCLImage],
    destination: Union[numpy.ndarray, pyclesperanto_prototype._tier0._pycl.OCLArray, pyopencl._cl.Image, pyclesperanto_prototype._tier0._pycl._OCLImage] = None,
    sigma_x: float = 0,
    sigma_y: float = 0,
    sigma_z: float = 0,
) -> Union[numpy.ndarray, pyclesperanto_prototype._tier0._pycl.OCLArray, pyopencl._cl.Image, pyclesperanto_prototype._tier0._pycl._OCLImage]
Docstring:
Computes the Gaussian blurred image of an image given sigma values
in X, Y and Z. 

Thus, the filter kernel can have non-isotropic shape.

The implementation is done separable. In case a sigma equals zero, the 
direction is not blurred. 

Parameters
----------
source : Image
destination : Image, optional
sigma_x : Number, optional
sigma_y : Number, optional
sigma_z : Number, optional

Returns
-------
destination

Examples
--------
>>> import pyclesperanto_prototype as cle
>>> cle.gaussian_blur(source, destination, sigma_x, sigma_y, sigma_z)

References
----------
.. [1] https://clij.github.io/clij2-docs/reference_gaussianBlur3D
File:      ~/Libraries/miniconda3/envs/proto/lib/python3.11/site-packages/pyclesperanto_prototype/_tier1/_gaussian_blur.py
Type:      function

gpu_image = cle.push(image)
gpu_blurred = cle.create_like(gpu_image)
cle.gaussian_blur(gpu_image, gpu_blurred, sigma_x=5, sigma_y=5)
blurred = cle.pull(gpu_blurred)

imshow(blurred)

/home/stephane/Libraries/miniconda3/envs/proto/lib/python3.11/site-packages/skimage/io/_plugins/matplotlib_plugin.py:149: UserWarning: Float image out of standard range; displaying image with stretched contrast.
  lo, hi, cmap = _get_display_range(image)

<matplotlib.image.AxesImage at 0x7f2b21949210>

Exercise 1: Apply a threshold Otsu

Now that we have a blurred image, we can apply an automatic thresholding method to segment the image. We will use the Otsu thresholding method for that. Push the blurred image into the GPU-memory, create an empty image as a destination output and apply the thresholding method threshold_otsu.

#TODO

imshow(binary)

<matplotlib.image.AxesImage at 0x7f2b2282bd50>

We have now done image processing on the GPU, Congratulations!

Advance syntax

What we saw so far is the basic syntax, which is pedagogically useful as well as meant for interoperability with other languages. However, for python it can be quite verbose and redondant, and does not take advantage of the python language features. The advance syntax is meant to be more pythonic, and should be more familiar to python users.

Let’s start again with our gaussian blur example

imshow(image)

<matplotlib.image.AxesImage at 0x7f2b21a5bd50>

Now lets apply the same operation (create, push, process, pull and display) on a 2D image and let’s see the difference!

gpu_blurred = cle.gaussian_blur(image, sigma_x=5, sigma_y=5)

print(type(gpu_blurred), gpu_blurred.shape, gpu_blurred.dtype)
gpu_blurred

<class 'pyclesperanto_prototype._tier0._pycl.OCLArray'> (254, 256) float32

cle._ image shape (254, 256) dtype float32 size 254.0 kB min 34.26861 max 237.56865

cle.gaussian_blur?

Signature:
cle.gaussian_blur(
    source: Union[numpy.ndarray, pyclesperanto_prototype._tier0._pycl.OCLArray, pyopencl._cl.Image, pyclesperanto_prototype._tier0._pycl._OCLImage],
    destination: Union[numpy.ndarray, pyclesperanto_prototype._tier0._pycl.OCLArray, pyopencl._cl.Image, pyclesperanto_prototype._tier0._pycl._OCLImage] = None,
    sigma_x: float = 0,
    sigma_y: float = 0,
    sigma_z: float = 0,
) -> Union[numpy.ndarray, pyclesperanto_prototype._tier0._pycl.OCLArray, pyopencl._cl.Image, pyclesperanto_prototype._tier0._pycl._OCLImage]
Docstring:
Computes the Gaussian blurred image of an image given sigma values
in X, Y and Z. 

Thus, the filter kernel can have non-isotropic shape.

The implementation is done separable. In case a sigma equals zero, the 
direction is not blurred. 

Parameters
----------
source : Image
destination : Image, optional
sigma_x : Number, optional
sigma_y : Number, optional
sigma_z : Number, optional

Returns
-------
destination

Examples
--------
>>> import pyclesperanto_prototype as cle
>>> cle.gaussian_blur(source, destination, sigma_x, sigma_y, sigma_z)

References
----------
.. [1] https://clij.github.io/clij2-docs/reference_gaussianBlur3D
File:      ~/Libraries/miniconda3/envs/proto/lib/python3.11/site-packages/pyclesperanto_prototype/_tier1/_gaussian_blur.py
Type:      function

If we look at the operation signature, we can observe that

• the destination output can be optional as it is None by default

• the source input can also be a numpy.ndarray, so not a GPU array

• the operation returns the destination output

This however does not change the process, if a numpy array is passed as input, it will be pushing it to the GPU, the operation will be performed. If no destination output is provided, we will automatically create a new GPU array and return it. Finally, the display still requires to pull the data back from the GPU in order to be readable.

TLDR: The librairy can consume numpy arrays directly and it will manage automatically the data transfer to and from the GPU for you (in most cases).

Array operators

Mathematics behing at the base of image processing, mathematical operators are also part of the image processing toolbox. A quick search of the available methods in the librairy reveals several mathematical operations:

• add_images

• binary_and

• add_image_and_scalar

• etc.

These operations can be use the same way as any other filters, providing an input image and retrieving an output image. It is also possible to use them as data operators -, +, *, /, etc. - on the image data. These operators follow the same usage as numpy arrays API.

print("minimum intensity:", gpu_blurred.min())
print("maximum intensity:", gpu_blurred.max())

minimum intensity: 34.26861
maximum intensity: 237.56865

gpu_blurred * gpu_binary

cle._ image shape (254, 256) dtype float32 size 254.0 kB min 0.0 max 237.56865

gpu_blurred - 50

cle._ image shape (254, 256) dtype float32 size 254.0 kB min -15.731392 max 187.56865

Exercise 2: Normalise the image gpu_blurred using operators

Inorm = (I - Imin) / (Imax - Imin)

#TODO

<class 'pyclesperanto_prototype._tier0._pycl.OCLArray'> (254, 256) float32

cle._ image shape (254, 256) dtype float32 size 254.0 kB min 0.0 max 1.0

Build processing pipeline

It is unlikely that you simply use one operation on your image. Usually, you want to combine multiple operations to achieve a certain goal. It is usually in this kind of cases that GPU processing really shines.
As for now, we always push and pull data to send it to the GPU and get it back. This is however not very efficient. Once the data is pushed to the GPU, you want to keep it there and perform as many operations as possible on it before having the need to pull it back.

Let’s see how we can do that for a basic pipeline:

gpu_image = cle.push(image)
blurred1 = cle.gaussian_blur(gpu_image, sigma_x=3, sigma_y=3)
blurred2 = cle.gaussian_blur(gpu_image, sigma_x=10, sigma_y=10)
difference = cle.subtract_images(blurred1, blurred2)
difference

Although we did three distinct operations, only one push/pull operation was performed. And using these three operations, we were able to perform a more complexe algorithm : the difference of gaussian filter.

Question : How many relative memory do I need to run this algorithm on my GPU ?

Tiers list of operations

Of course, we are not going to call the three operations everytime we want to apply the difference_of_gaussian filter. We already have this operations already implemented in the pyclesperanto_prototype library.

You can have a look at it’s implementation here as well as see how to use it with the ? operator.

cle.difference_of_gaussian?

You can have a complete list of the available operations in the documentation, or by browsing the source code of the library. You should also take advantage of the auto-completion feature of your IDE to discover the available operations. As the librairy grows we, little by little, add more and more advance and complexe workflows simplifying some redondant processing.

def list_operations(search_term = None):
    ops = cle.operations(search_term)
    for name in ops:
        op = ops[name]
        func = cle.operation(name)

        if hasattr(func, 'fullargspec'):
            print(name + "(" + str(func.fullargspec.args).replace('[','').replace(']','').replace('\'','') + ")")
        else:
            print(name)

list_operations()