Basic image manipulation

Previously, you've learned to load and save medical image with MedPy. That was really easy, wasn't it? But wait until you'll see what power you've just obtained. Take a closer look at the following examples of basic image manipulation.

Thresholding

Thresholding is an often used image manipulation task. The goal is to create a binary mask that selects certain voxels based on their intensity values. To e.g. mark all voxels with intensity values greater than 100, you can simply call

mask = image > 100

This results in a binary image of the same shape as the original image with True-values where the original image's intensity is greater than 100 and False values everywhere else.

image.shape
> (91, 109, 91)
mask.shape
> (91, 109, 91)
mask.dtype.type
> numpy.bool_

Assuming you have an upper as well as a lower threshold, you can use

mask = (image > 100) & (image < 200)

Note the brackets around the expressions, which are required since the logical-& is stronger than the >.

To obtain a mask for all voxels with a certain intensity value, e.g. 0 in this case, call

mask = 0 == image

I use here the constant-first notation, but

mask = image == 0

would be equivalent, only slightly less readable and fail-safe.

I've talked about the power of numpy and we are actually already using it all the time. Now we will inclue the module itself and use one of its convenient functions to count the number of True values in the mask

import numpy
numpy.count_nonzero(mask)
> 60

Working with a mask

background_mask = image < 20
background_mask.dtype
> numpy.bool
image[background_mask] = 0
image[background_mask].shape
> (364739,)

Region growing

Let us implement a numpy-style n-dimensional region growing algorithm.

import numpy
from scipy.ndimage import label
def region_growing(img, seed, minthr, maxthr, structure = None):
	img[seed] = minthr
	thrimg = (img < maxthr) & (img >= minthr)
	lmap, _ = label(thrimg, structure = structure)
	lids = numpy.unique(lmap[seed])
	region = numpy.zeros(img.shape, numpy.bool)
	for lid in lids:
		region |= lmap == lid
	return region

Neat, eh? Lets get through it step by step. First we define the function's signature.

def region_growing(img, seed, minthr, maxthr, structure = None):

It expects an image that is a numpy ndarray, a seed binary image of the same shape that contains a True-value at each seed point, a minimal threshold value, a maximum threshold value and, optionally, a structure element that describes the connectedness / neighbourhood.

img[seed] = minthr

First we set all seed points in the image to the minthr value. We will go into this later.

thrimg = (img < maxthr) & (img >= minthr)

We threshold the image and obtain a binary image with True-values for all voxels that are larger or equal than minthr and smaller than maxthr. These voxels are candidates for our segmentation.

lmap, _ = label(thrimg, structure = structure)

Now we use the labelling from scipy.ndimage, which labels all unconnected (where connectedness is defined by our structure element) binary objects in the threshold image with different label ids. The result is hence an unsigned integer image.

lids = numpy.unique(lmap[seed])

Now we mask the label map image with our seed binary mask and obtain a list of the label ids at the position of the seeds. Of course, it would be problematic, if the seeds would lie in the background area. That's where our first line comes into play: We assigned the minthr value to these voxels, such that they are guaranteed to have been segmented in the segmentation step.

region = numpy.zeros(img.shape, numpy.bool)
for lid in lids:
	region |= lmap == lid

We now have the ids of all labels in which lies at least one of our seed points. We simply extract their locations from the label map and combine them into a segmentation region image by the logical-or operator.

return region

Finally, we return the obtained segmentation region.

Done. Note the way we use a combination of masks, labels and existing function for a fast, robust and n-dimensional region growing implementation without ever attempting any iteration over the image, but rather implicitly making use of the numpy broadcasting rules. This is very important if your goal is speed as often the case in the processing of large images.