Slice by Slice Adaptive Histogram Equalization

Overview

Most SimpleITK filters can only operate on 2 or 3 dimensional images, with the exception of filters such as ExtractImageFilter, PasteImageFilter, SliceImageFilter and JoinSeriesImageFilter. However, SimpleITK (by default) supports upto 5 dimensional images. A high dimensional image can be processed by extracting 2 or 3 dimensional images, then either using the JoinSeriesImageFilter to join the sub-volumes together or the PasteImageFilter to copy the results back to the original image. Possible reasons include when the z direction spacing is too great, or for computation or memory efficient reasons. Additionally, it may be desired to process a volume (3d or higher) as a sequence of 2 dimensional images.

In this example, the AdaptiveHistogramEqualizationImageFilter is used to processes a higher dimensional image as a sequence of two dimensional images.. When the filter is run only a single X, Y cross-section of the volume.

Both the Python and the C++ examples demonstrate a reusable “decorator” to wrap the SimpleITK AdaptiveHistogramEqualization procedure to process the input image by 2d slices. A function decorator is a function which takes a function as an argument and returns a modified or wrapped function. The decorators are written generically, so they can work with other SimpleITK procedures or a custom procedure. The decorators wrap the procedure to accept a 2, 3, 4 or 5 dimensional image as input. The pasting approach is used to provide a memory efficient approach.

The process of extracting a slice, processing, and then pasting back to the original image is straight forward. However, to create a reusable decorator requires advanced language specific features. Additionally, to efficiently do pasting in place execution is done in C++, and sliced indexed assignment is used in Python.

Code

PythonC++

#!/usr/bin/env python

""" Example demonstrating a slice-by-slice function decorator """

import sys
import itertools
from functools import wraps
import SimpleITK as sitk


def slice_by_slice_decorator(func):
    """
    A function decorator which executes func on each 3D sub-volume and
    *in-place* pastes the results into the input image. The input image type
    and the output image type are required to be the same type.

    :param func: A function which take a SimpleITK Image as it's first
    argument and returns an Image as results.

    :return: A decorated function.
    """

    iter_dim = 2

    @wraps(func)
    def slice_by_slice(image, *args, **kwargs):

        dim = image.GetDimension()

        if dim <= iter_dim:
            #
            image = func(image, *args, **kwargs)
            return image

        extract_size = list(image.GetSize())
        extract_size[iter_dim:] = itertools.repeat(0, dim - iter_dim)

        extract_index = [0] * dim
        paste_idx = [slice(None, None)] * dim

        extractor = sitk.ExtractImageFilter()
        extractor.SetSize(extract_size)

        for high_idx in itertools.product(
            *[range(s) for s in image.GetSize()[iter_dim:]]
        ):

            # The lower 2 elements of extract_index are always 0.
            # The remaining indices are iterated through all indexes.
            extract_index[iter_dim:] = high_idx
            extractor.SetIndex(extract_index)

            # Sliced based indexing for setting image values internally uses
            # the PasteImageFilter executed "in place".  The lower 2 elements
            # are equivalent to ":". For a less general case the assignment
            # could be written as image[:,:,z] = ...
            paste_idx[iter_dim:] = high_idx
            image[paste_idx] = func(extractor.Execute(image), *args, **kwargs)

        return image

    return slice_by_slice


if len(sys.argv) < 3:
    print("Usage: SubDimensionProcess inputImage outputImage", file=sys.stderr)
    sys.exit(1)

inputImage = sitk.ReadImage(sys.argv[1])

# Decorate the function
adaptive_histogram_equalization_2d = slice_by_slice_decorator(
    sitk.AdaptiveHistogramEqualization
)

adaptive_histogram_equalization_2d(inputImage, radius=[20] * 2, alpha=0.3, beta=0.3)

sitk.WriteImage(inputImage, sys.argv[2])