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¶
#!/usr/bin/env python
from __future__ import print_function
import SimpleITK as sitk
import sys
import itertools
from functools import wraps
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 "inplace".
# 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])
#include "SimpleITK.h"
#include "sitkExtractImageFilter.h"
#include <cstdlib>
#include <iostream>
namespace sitk = itk::simple;
using itk::simple::operator<<;
// Forward declaration to specialize the implementation with the function's
// argument and returns types.
template <class T> struct SliceBySliceDecorator;
/* \brief A function decorator to adapt an function to process an image as a
* sequence of 2D (slices or) images.
*
* For ease of use the makeSliceBySlice procedure should be used to construct
* the SliceBySliceDecorator.
*
* The overloaded function call operator(), enable objects to be used as
* functions.
*
* The return image is the first argument modified with the slice by slice
* results of the f function.
*/
template <class R, class ImageArg, class... Args>
struct SliceBySliceDecorator<R(ImageArg, Args...)>
{
using FunctionType = std::function<R(ImageArg, Args...)>;
explicit SliceBySliceDecorator(FunctionType f) : f_(std::move(f)) {}
R operator()(sitk::Image &image, Args... args)
{
const auto image_size = image.GetSize();
unsigned int dim = image.GetDimension();
if (dim <= iter_dim)
{
// If no sub-dimension extraction is required then directly run the
// function on the input and replace it.
image = f_(image, args...);
return image;
}
std::vector<unsigned int> extract_size = image.GetSize();
std::fill(extract_size.begin() + iter_dim, extract_size.end(), 0);
std::vector<int> extract_index(dim, 0);
// The extract filter is used to extract a 2D image from the higher
// dimensional input image.
sitk::ExtractImageFilter extractor;
extractor.SetSize(extract_size);
// The paste filter places the processed slice back into the original image
// at the same location where the extraction occurred. The default
// value of DestinationSkipAxes is [ false, false, true, ... ], which is
// correct for the situation of preserving the first dimensions and
// collapsing the remainder.
sitk::PasteImageFilter paster;
paster.SetSourceSize(std::vector<unsigned int>(
extract_size.begin(), extract_size.begin() + iter_dim));
while (static_cast<unsigned int>(extract_index.back()) < image.GetSize().back())
{
extractor.SetIndex(extract_index);
// Store the results of the function as a r-value, so that the
// paste filter will run "in place" and reuse the buffer for output.
sitk::Image &&temp_image = f_(extractor.Execute(image), args...);
paster.SetDestinationIndex(extract_index);
image = paster.Execute(image, temp_image);
// increment extraction index
++extract_index[iter_dim];
for (unsigned int e = iter_dim; e + 1 < dim; ++e)
{
// if the index element is beyond the size, propagate to next element.
if (static_cast<unsigned int>(extract_index[e]) > image.GetSize()[e])
{
extract_index[e] = 0;
++extract_index[e + 1];
}
}
}
return image;
}
FunctionType f_;
constexpr static unsigned int iter_dim = 2;
};
/** Construct a function object to wrap a function f to process a volume slice by
* slice with the f function.
*
* The f function must take an sitk::Image object as the first parameter and
* return an sitk::Image.
*/
template <class R, class... Args>
SliceBySliceDecorator<R(Args...)> makeSliceBySlice(R (*f)( Args...))
{
using DecoratorType = SliceBySliceDecorator<R(Args...)>;
return DecoratorType(typename DecoratorType::FunctionType(f));
}
int main( int argc, char *argv[])
{
if ( argc < 3 )
{
std::cerr << "Usage: " << argv[0] << " <inputImage> <outputImage>" << std::endl;
return 1;
}
sitk::Image tempImage = sitk::ReadImage(argv[1]);
// The parameters to the filter are hard coded to simplify the example.
float alpha = 0.3f;
float beta = 0.3f;
std::vector<unsigned int> radius(2, 20);
// The return type of the function decorator is complex, so the auto type is
// used for brevity.
auto AdaptiveHistogramEqualization2D = makeSliceBySlice(sitk::AdaptiveHistogramEqualization);
// Execute the decorated function.
tempImage = AdaptiveHistogramEqualization2D(tempImage, radius, alpha, beta);
sitk::WriteImage(tempImage, argv[2]);
return 0;
}