Image Segmentation, Registration and Characterization in R with SimpleITK

Many types of medical and scientiﬁc experiments acquire raw data in the form of images. Various forms of image processing and image analysis are used to transform the raw image data into quantitative measures that are the basis of subsequent statistical analysis. Inthis article we describe the SimpleITK R package. SimpleITK is a simpliﬁed interface to the insight segmentation and registration toolkit ( ITK ). ITK is an open source C++ toolkit that has been actively developed over the past 18 years and is widely used by the medical image analysis community. SimpleITK provides packages for many interpreter environments, including R . Currently, it includes several hundred classes for image analysis including a wide range of image input and output, ﬁltering operations, and higher level components for segmentation and registration. Using SimpleITK , development of complex combinations of image and statistical analysis procedures is feasible. This article includes several examples of computational image analysis tasks implemented using SimpleITK , including spherical marker localization, multi-modal image registration, segmentation evaluation, and cell image analysis.


Introduction
Images are an important source of quantitative information in many fields of research and many industrial, medical and scientific applications. The processing required to transform pixel or voxel data into quantitative information used in subsequent statistical analysis is as varied as the imaging instruments used to acquire image data, the scenes being imaged and the form of information required. Instruments in typical use include basic digital cameras in scenarios like security and wildlife monitoring, digital cameras associated with an ever increasing variety of microscopes, and medical imaging equipment such as magnetic resonance imaging (MRI), computed tomography (CT), ultrasound and X-ray. Other instruments, such as LIDAR, produce gridded data that may be useful to regard as an image in some scenarios. This small list of examples serves to illustrate the range of data types typically acquired, before considering the range of scenes being imaged. Typical digital cameras are capable of acquiring two-dimensional scenes, possibly with color channels and possibly with time series. MRI is capable of acquiring a range of three-dimensional (3D) time series and other, more specialized, forms including directional information in 3D. Microscopes are capable of acquiring time series of 3D images with multiple fluorescent channels, and such microscopes are being coupled to automated sample preparation systems to produce automated screening platforms capable of generating massive amounts of images from complex experiments.
The type of information extracted from an image may range from a detection (e.g., presence of an intruder or an animal) to a count (e.g., number of cells or brain lesions) to more complex characterizations of objects (e.g., measures of length/area/volume, spatial distributions, brightness) to tracks of moving objects in time series. There are numerous options, and the difficulty varies considerably with scene complexity and consistency.
Creating a new image analysis process typically requires experimentation with a variety of approaches, each of which will combine several computational steps or components. It is a considerable advantage to have on hand an image analysis framework with a large number of functional components capable of operating on multi-dimensional data when embarking on such an endeavor. The Insight Segmentation and Registration Toolkit (ITK; Johnson, McCormick, Ibáñez, and The Insight Software Consortium 2013) is such a framework, offering thousands of components for image input and output, image filtering, image segmentation and image registration. Originally designed for medical images, ITK has been successfully used in varied domains with data ranging from cell images to remote sensing images.
In this article we introduce the SimpleITK package (The Insight Software Consortium 2018) for R(R Core Team 2018), which allows development of image analysis pipelines using R and a simplified interface to ITK.
ITK has two features distinguishing it from most image analysis toolkits: • Complex operations are accomplished by using filter pipelines.
Filter pipelines are the architecture used in ITK to accomplish complex tasks, such as segmentation, via lazy evaluation. A pipeline of filters is constructed to implement a series of computational steps, however no computation is performed until a call is made to the Update method of the last filter. The request is propagated back along the pipeline and only the required processing is performed. Pipelines also facilitate processing of datasets that are too large for core memory via streaming.
• An image occupies a region in physical space and operations should take this into account.
The center of every voxel in an ITK image is located at a specific point in space with respect to a coordinate system, the axes of which are not necessarily orthogonal. Distances between adjacent voxel centers are in physical units and are not necessarily isotropic.
ITK offers a high degree of customizability which, unfortunately, leads to a steep learning curve for developers and requires significant C++ experience.

SimpleITK
By making several judicious design decisions, SimpleITK is able to provide a simplified interface to ITK while maintaining most of the functionality available in the original toolkit. The key design decisions were based on the following observations obtained via surveys: most users analyze two-dimensional (2D), three-dimensional (3D), and possibly four-dimensional images (3D plus time); most users do not take advantage of the features of the pipeline approach; non-C++ developers prefer a procedural interface over an object oriented one.
As a consequence, SimpleITK was designed to support 2D and 3D images, with an optional configuration flag for 4D support. The pipeline architecture is hidden in favor of immediate execution, and both, an object oriented and a procedural interface, have been implemented. Note that while the image dimensions are restricted, voxels can be vectors of arbitrary length containing a multi-dimensional measure, such as color, fiber orientation or spectra.
Just like ITK, SimpleITK is implemented in C++. Unlike ITK, most of the code is generated automatically using a JavaScript Object Notation (JSON) template structure. A detailed description of the design and automated code generation process is given in Lowekamp, Chen, Ibáñez, and Blezek (2013). The design of SimpleITK also allows for easy automated wrapping of the C++ implementation for several interpreted languages including Python (Rossum et al. 2011), Perl (Christiansen, Foy, Orwant, and Wall 2012), C#, Lua (Ierusalimschy 2016), Ruby (Thomas, Fowler, and Hunt 2009), Tcl (Tcl Core Team 2017), and R, the focus of this article.
Automated generation of interface code is critical for long term viability of SimpleITK due to its dependence on ITK, which is very large and continually evolving.

Licensing
SimpleITK is distributed under the Apache 2.0 license. The full license is available at http: //www.simpleitk.org/SimpleITK/project/license.html.

The SimpleITK R package
Multiple layers of automated code generation and dependencies on external tools make building SimpleITK a challenge and cause problems distributing packages via traditional mechanisms, such as at the Comprehensive R Archive Network (CRAN). Two approaches for Linux and OSX systems are described below.

Documentation
Documentation is automatically converted from the doxygen extracted C++ class documentation. It offers a non-standard starting point for the R developer. Most important details are to be found in the help for the class interface, rather than the help for the procedural interface. For example ?Cast will display only the most basic information concerning usage (useful for determining argument names and order) while details of functionality are available in ?CastImageFilter. This division is common across SimpleITK and is shared with the C++ documentation. The C++ documentation is structured to describe the classes and associated methods, and thus does not fit into the R function documentation design. The current approach maps class methods to function arguments. The SimpleITK C++ documentation (The Insight Software Consortium 2016b) is the canonical source of information concerning available methods and classes.
Beyond the API documentation described above the toolkit also maintains general, language agnostic, documentation on read-the-docs https://simpleitk.readthedocs.io/. This documentation covers installation, fundamental concepts specific to SimpleITK's image and transformation elements, common API conventions, frequently asked questions and short example programs. Additional resources include a Jupyter notebook repository illustrating complete image-analysis workflows in Python and R (Yaniv, Lowekamp, Johnson, and Beare 2018, ; https://github.com/InsightSoftwareConsortium/SimpleITK-Notebooks) and a discourse discussion forum for users to post questions (https://discourse.itk.org/).

Building from source
Greater control of the build process, as well as the ability to generate other wrapper packages and run testing procedures, is available with a source build.
Building of the R package involves two code generation steps. In the first step the C++ classes are generated using scripts written in the Lua programming language. These classes are described by a combination of template and JSON files describing each of the classes. Some additional classes are implemented directly in C++. Once all of the C++ source is created, SWIG is used to generate the R bindings for all classes.
The current build process requires CMake and Git to fetch and build dependencies. The additional components on which the build process depends include SWIG, Lua, and ITK. For a fully automated build process without any pre-installed components one selects the project's SuperBuild CmakeLists file without making any change to the CMake settings. The build process is documented at The Insight Software Consortium (2016a). In brief, the process is: • Clone the source code: $ git clone http://itk.org/SimpleITK.git • Create a build folder and run cmake on the SuperBuild subdirectory: $ mkdir SimpleITK-build $ cd SimpleITK-build $ cmake ../SimpleITK/SuperBuild $ make The build process fetches and builds all dependencies. A successful build produces a binary R package in SimpleITK-build/Wrapping/R that can be distributed and installed.

Wrapper generation with SWIG
The software interface generator, SWIG (Beazley et al. 1996), was originally designed to generate Python interfaces to C code, but now supports many interpreted environments, including R. This makes SWIG a natural choice as the interface generator for SimpleITK, which aims to support multiple interpreted programming languages. However SWIG's popularity in the R world is probably lower than Rcpp (Eddelbuettel and François 2011), which provides an API layer for C++/R rather than an interface generator. The SWIG support for R has been extended to support SimpleITK development, specifically support for C++ stl vectors (automatic conversion to and from R vectors). These developments potentially make SWIG an interesting alternative to Rcpp for projects targeting multiple programming languages in addition to R.
Objects in SimpleITK, such as images and filters, are external references to C++ objects. SWIG generated binding code uses a combination of S4 classes and a reference class style structure to wrap these objects which are then managed via the R garbage collection system.

Multi-threaded components
There are a number of architectural features in ITK to make construction of multi-threaded filters simpler. As a result a substantial proportion of SimpleITK components are able to take advantage of multicore processors and multiple processor configurations. This includes all operations that operate on single pixels, such as image arithmetic, many kernel filters and parts of the registration framework. By default the number of threads a component uses is equal to the number of cores available. The number of threads used by a class can be queried and set as shown in the following code snippet.

Image files, image objects and meta-data
SimpleITK can read and write many image file formats, including traditional 2D formats such as jpeg, bmp, png and tiff, as well as medically oriented, 3D formats including Digital Imaging and Communications in Medicine (DICOM), analyze, nifti, nrrd, metaimage and others. The image format is automatically determined during reading and writing.
Images in SimpleITK are external references, that is pointers to C++ objects. They are created by loading data from files, calls to constructor methods or conversion from R arrays. Image meta-data, such as voxel type, voxel spacing, and image orientation may be queried and set using method calls, as follows: R> img <-ReadImage("phantom.dcm") R> img$GetSize() This DICOM image is a single 512×512 slice, with voxel dimensions 0.65625×0.65625×1mm.
The presence of slice thickness information in the DICOM file leads to creation of a 3D image, even though there is only a single slice. Image formats, such as png, which typically do not have thickness information will be interpreted as 2D images.

Image extent and physical coordinates
Relationships between image data and the patient are critical when working with medical images. It is usually important to know the precise size of a voxel, and sometimes exactly where it was located in the scanner. Such knowledge is clearly needed when measuring sizes of pathology, such as tumors, that are expected to change size over time. It is also critical in applications where images are combined, to ensure that only legal operations are performed -for example, performing masking with images of the same size but different orientations is likely to cause errors. SimpleITK provides image orientation and image voxel spacing information that enable these operations and ensure that illegal operations are detected.
The key components are image origin, voxel spacing and axis orientation. The following code snippet illustrates the creation of a 2D axis aligned image with non-zero origin and nonisotropic pixel spacing and demonstrates mapping between voxel and physical units. Looking at the returned indexes we see that the pixel extent in physical space indeed starts half a pixel from its center and that the indexes in SimpleITK are zero-based.
SimpleITK has no inherent image display capabilities. A convenience Show function is provided, displaying images by writing them to a temporary location in nifti format and then invoking an external viewer. By default this function will attempt to invoke the Fiji/ImageJ viewer. Invoking other viewers is facilitated by setting SITK_SHOW_COMMAND as environment variable to indicate a specific viewer. Finer grained control is possible, allowing one to specify a viewer for color images and for 3D images by setting the SITK_SHOW_COLOR_COMMAND and SITK_SHOW_3D_COMMAND environment variables.

Conversion to and from arrays
Images can be converted to native R arrays using as.array and arrays can be converted to images using as.image. Image meta-data can be copied from one image to another or set manually. Note that as.array includes a drop option to discard unit dimensions.

Indexing and slicing
Unlike R indexing, SimpleITK functions use zero-based indexing, so you will need to account for this when working with SimpleITK images. The exception to the rule is slicing. This is a native R operation and SimpleITK supports it using R's one-based indexing: Slicing allows complex manipulation of images using standard syntax, returning new images. Cropping, flipping, etc. is easy to achieve using slicing operations: • Extracting one corner: R> img <-Image(10, 10, "sitkUInt8") R> imgCorner <-img[1:5, 1:5] R> imgCorner$GetSize() [1] 5 5 • Simple reflection: • Accessing a single pixel: • Illegal operation -remove central columns producing non-uniform spacing: Error in img [-c(3:5), ]: X spacing is not uniform Image slicing is designed to preserve image constraints, and is thus slightly less flexible than standard R array indexing. A slicing operation must produce images with uniform voxel spacing along each image dimension. The concept of empty arrays, produced by indexing operations like ar[0, 1:2], does not translate directly to image objects, and an error is raised. A voxel is converted to a native R object when a slicing operation produces a single voxel result (rather than returning a single voxel image). Slicing of multi-component images is possible, with multi-component voxels being returned as R vectors. Array-based indexing is not yet supported.

Error messages
The reference class style of method call, combined with the need for method overloading, leads in some cases to unhelpful error messages, such as when invoking a non-existent method: R> img <-Image(10, 10, "sitkUInt8") R> img$SetVoxel(c(5, 5), 255)

Error in callNextMethod() : node stack overflow
Or passing combinations of arguments that the method overloading code is unable to resolve: R> th <-Threshold(img, "pixtype") Error in f(...): could not find function "f" In the first example we should be using SetPixel while in the latter case there is no function call version of Threshold with image and character arguments.
Finally, configuration of the external viewer often causes problems the first time one uses SimpleITK. The primary reason being that the default viewers, ImageJ/Fiji, are not installed and the environment variable specifying an alternative viewer was not set as described above.

SimpleITK computational components
SimpleITK reduces the learning curve compared to ITK by simplifying many of the programming aspects. Unfortunately, mastering all of the available functionality still requires an effort, primarily due to the size of the computational framework.
It is often difficult to identify components of interest due to differences in naming conventions between toolkits. ITK uses a module hierarchy for code that can help make the range of components more comprehensible. Not all of the classes are directly accessible from SimpleITK, but most are or will eventually be. Infrastructure modules have been left out for clarity. The key modules and their contents are described below.

Filtering
This is the largest module and contains numerous classes that modify voxel values in some way. It contains the following categories, in alphabetical order: • AnisotropicSmoothing: Gradient and curvature anisotropic diffusion for scalar and vector images.
• AntiAlias: Reduce aliasing artifacts when displaying binary volume.
• BiasCorrection: The N4 algorithm for bias field inhomogeneity in MR images.
• BinaryMathematicalMorphology: Specialized mathematical operations for binary images.
• Colormap: Colormaps for creating overlays and display.
• DisplacementField: Tools for processing displacement images produced by registration.
• DistanceMap: Signed and unsigned distance transforms using various algorithms.
• FastMarching: Classes to compute geodesic distances on images using the fast marching algorithm.
• FFT: A range of fast Fourier transform classes.
• ImageLabel: Manipulation of label images (outputs of connected component analysis).
• LabelMap: Manipulation of run-length-encoded objects from connected component analysis.
• Thresholding: Threshold estimation using a range of histogram-based approaches.

Segmentation
Image segmentation filters produce output images in which voxel values indicate class membership. The relevant ITK modules are: • Classifiers: Bayesian, K-means voxel classifiers.
• ConnectedComponents: Label objects in a binary image (aka particle analysis).
• DeformableMesh: Mesh-based segmentation in which mesh deformation is driven by image forces.
• LabelVoting: Various schemes for combining label images, including STAPLE, voting. Hole filling with voting.
• LevelSets: An extensive framework for segmentation using the level set methodology. Includes geodesic active contours, shape priors and others.
• RegionGrowing: Various combined threshold and connectivity approaches to segmentation.
• Watersheds: Several algorithms for watershed transforms, including marker-based options.

Registration
Image registration is the process of estimating a spatial transformation that makes two images similar according to some measure. It is usually structured as an optimization problem. There are numerous choices available for similarity measures, optimizers, and transformation functions. For intensity-based registration, SimpleITK makes many of these available via a single framework class, with options for callback functions to track optimizer progress. The available choices include: • Similarity metrics: Correlation, mutual information (Mattes and joint histogram), Demons, mean squares, and the ANTS neighborhood correlation.
Classes implementing transforms, optimizers and interpolators are in core modules, rather than the registration module, as they are used in other scenarios.
Additional algorithm specific implementations are also available. These include the implementation LandmarkBasedTransformInitializer, estimating the transformation which minimizes the distance between two point sets with known correspondences (transformation can be rigid, affine or B-spline); and several Demons-based intensity-based algorithms for estimating the dense deformation field between images, Demons, DiffeomorphicDemons, SymmetricForcesDemons and FastSymmetricForcesDemons.

Case studies
The addition of SimpleITK to the R development environment enables the rapid development of end-to-end image characterization and statistical analysis tools. This combination benefits from the large choice of image operators available in SimpleITK and the extensive range of existing R packages for feature extraction and statistical analysis. In addition, the use of SimpleITK in an interpreted environment offers a reduction in development cycle time by removing the compile stage from the usual change-compile-test cycle. In this section we will illustrate SimpleITK via several example case studies, with the aim of providing an introduction to both, the slightly unusual syntax of a SWIG interface and some of the extensive capabilities of SimpleITK. These examples illustrate only a small proportion of the available classes but provide an introduction to several components that are useful in many different scenarios.

Spherical marker localization
Alignment between a patient and volumetric images of that patient is critical in many forms of surgical/medical intervention. Spherical markers that are visible on the patient and in the volumetric images are frequently used to aid registration in computer assisted intervention. Before the markers can be used for alignment they need to be localized in the image coordinate system as described in Yaniv (2009). There are a variety of approaches for spherical marker localization. In our case we will perform edge detection using SimpleITK, and use R to obtain a least-squares estimate of the sphere's location and radius.
Our input image was acquired using a cone-beam CT system and contains a single sphere. The image is non-isotropic, a 40 × 25 × 11 volume with a 0.44 × 0.44 × 0.89mm spacing, and has a high dynamic range of intensity values, [−32767, −25481]. Figure 1 shows the volume containing the metallic sphere and the result of performing edge detection on it.
The localization code, below, employs edge detection with parameters selected to match the anisotropic image voxel spacing.
First, load the image.
R> library("SimpleITK") R> sphere_image <-ReadImage("sphere.mha") Second, perform 3D edge detection to produce a binary edge map using a Canny edge detector. Third, convert the edge image to an array and determine the physical coordinates of each edge voxel using the built-in image method TransformIndexToPhysicalPoint. Note that indexes are adjusted to zero-based to match C++ standards.

Intensity-based image registration
Registration is the process of computing the transformation that relates corresponding points in two different coordinate systems. It is a key component in many medical and non-medical applications, with a large number of algorithms described in the literature (Oliveira and Tavares 2014;Zitová and Flusser 2003).
Intensity-based image registration estimates the transformation that aligns images to each other based on their intensity values. The task is formulated as an optimization problem with the optimized function indicating how similar the two images are to each other after applying a given spatial transformation. As an image consists of a discrete set of intensity values at grid locations and the transformation is over a continuous domain, this formulation also requires the use of interpolation. The terminology which we use here refers to one image as the fixed image and the other, which is being transformed to match the fixed image, as the moving image. The transformation whose parameters we are optimizing maps points from the fixed image coordinate system to the moving image coordinate system. Registration is discussed in more detail in Fitzpatrick, Hill, and Maurer Jr (2000), Yaniv (2008), or Goshtasby (2005).
The four components that one needs to specify in order to configure the registration framework in SimpleITK are: 1. Transformation type -global domain (linear) transformations such as AffineTransform or Euler3DTransform are available, or local domain (nonlinear) transformations such as BSplineTransform and DisplacementFieldTransform. This choice defines the set of parameters whose values are optimized.
2. Similarity function -a model of the relationship between the intensity values of the two images such as Correlation for affine relationship and MattesMutualInformation for a stochastic relationship. The function is expected to attain a local minimum when the images are correctly aligned. 3. Interpolator -method to use for interpolating intensity values at non-grid locations, such as sitkNearestNeighbor, sitkLinear or sitkHammingWindowedSinc. This choice often reflects a compromise between accuracy and computational complexity with the most common choice being sitkLinear.
4. Optimizer -algorithm used to reach the optimum of the similarity function. These range from simple evaluation of the similarity using a discrete grid in the parameter space, Exhaustive, to a limited memory method with simple constraints, L-BFGS-B. There can be a complex relationship between an optimizer and parameters of a transform, especially when transform parameters have very different units (e.g., radians and millimeters).
An optional additional component, an observer function, can be added to report registration progress. One can add multiple functions to the same observed event or different functions for each observed event. Observer functions written in R are readily utilized.
In the following example we align a patient's cranial CT to their T1 MRI scan. This data is available online and is the training data set provided as part of the Retrospective Image Registration Evaluation project (Fitzpatrick 2016 Figure 2 shows five axial slices from the data at different phases of the registration process.
The specific component selections for the task at hand are as follows. As both images were obtained from the same patient, the transformation between them is expected to be rigid. We use a Euler angle-based parameterization of the transformation. As the intensities of the two modalities are related via a stochastic relationship we use mutual information as the similarity function. To reduce the computational burden of computing the similarity function we use 1% of the voxels, selected via simple random sampling, leading to slight differences in final transform between runs. The other available option is to obtain them using a regular grid overlaid onto the image. We use linear interpolation to obtain intensity values at non-grid locations. Finally, we use the basic gradient descent optimizer as our optimization method.

Create the registration object and attach observers:
R> reg <-ImageRegistrationMethod() R> reg$AddCommand("sitkStartEvent", start_plot) R> reg$AddCommand("sitkIterationEvent", function() plot_values(reg)) Configure the registration object and execute: Two features of the registration framework that are specific to ITK and SimpleITK are the use of a so-called centered transform, CenteredTransformInitializer, and the automated estimation of parameter scales, SetOptimizerScalesFromPhysicalShift, for gradient-based optimizers. The centered transform performs rotation about the center of the fixed image, rather than the origin of the coordinate frame. The automated estimation of parameter scales deals with the complex relationship between transform terms with different units. For a comprehensive overview of the ITK registration framework we refer the interested reader to Avants, Tustison, Stauffer, Song, Wu, and Gee (2014).
Finally, even when the registration appears to have been performed successfully as validated by visual inspection (i.e., Figure 2) we always check the reason for the optimizer's termination. In our case we see below that we are likely dealing with premature termination, as the optimizer terminated because it reached the maximal number of iterations which we set. However the rate of improvement in similarity measure, illustrated in Figure 3 and based on information stored by the observer functions, decreased markedly after iteration 50. The solution is therefore unlikely to improve much with further iterations.

R> reg$GetOptimizerStopConditionDescription()
[1] "GradientDescentOptimizerv4Template:" [2] "Maximum number of iterations (100) exceeded." Rigid body registration, as illustrated in this example, is useful for aligning images of the same patient taken at similar times. This includes different modalities acquired during a study, such as different MR weightings, CT or PET as well as time series images, such as fMRI, where rigid body registration can be used to correct for patient movement.

Segmentation evaluation
Evaluation of segmentation algorithms applied to natural, medical, and biomedical images is most often done by comparing the output from the algorithm to those of human raters. It is common practice to derive an underlying reference segmentation by combining the annotations obtained from multiple raters. All raters can then be compared to the reference, which is useful when one rater is not known to be better than all others. Creating a combined reference is not straightforward, i.e., majority vote, when the raters are lay people such as when using crowd sourcing platforms as Amazon's mechanical Turk (Ipeirotis, Provost, and Wang 2010). Surprisingly, this is also not straightforward when the raters are domain experts such as radiologists interpreting medical images (Warfield, Zou, and Wells III 2004).
In this example we illustrate how to use the simultaneous truth and performance level estimation (STAPLE) algorithm (Warfield et al. 2004) to derive a reference segmentation. We then use this derived segmentation to evaluate the quality of the segmentation provided by each of our raters. As no single quality measure reflects the quality of the segmentation we generate a set of quality measures relating to overlap, over-and under-segmentation (false Figure 4: Deriving a reference segmentation from multiple raters using the STAPLE algorithm: (top row) manual segmentations performed by three radiologists; (bottom row) two additional segmentations derived from the expert segmentations to illustrate the effects of over-segmentation and segmentation with outliers. Last image is the derived reference segmentation obtained by the STAPLE algorithm.
positive, false negative), and distances between the boundaries of the segmented objects. The overlap and error scores between regions S and T are defined as follows: • Dice: 2 |S ∩ T | |S| + |T | .
• False positive: Scores for boundary distances are based on summaries (max, mean, median and standard deviation) of the distance between the contour defined by S and the closest point on the contour of T . These measures characterize difference in outline. When comparing two large regions, a large maximum distance between boundaries may result in a very small difference in overlap scores, for example if caused by a single isolated voxel. Whether this distinction is important is application dependent.
In our case we use a single slice from an abdominal CT scan in which three radiologists segmented a liver tumor. We added two additional segmentations with intentional errors resulting in over-segmentation and segmentation which contains an outlying region. Figure 4 visually illustrates the inputs to the STAPLE algorithm and the derived reference segmentation.
We next look at the code used to generate the reference segmentation and how it is used to evaluate segmentations. We start by loading the segmentations provided by our raters.
R> foreground_value <-1 R> threshold <-0.95 R> reference_segmentation_STAPLE_probabilities <-STAPLE(segmentations, + foreground_value) R> STAPLE_reference <-reference_segmentation_STAPLE_probabilities > + threshold Using the derived reference segmentation we can compare how each of the raters agrees with our reference segmentation. Note that in practice you would compare a new rater or a segmentation algorithm's performance to the derived reference segmentation. Two common ways to perform this evaluation include computation of overlap and boundary distance scores.
Computing the overlap scores is straightforward, simply provide the two segmentations as input to the LabelOverlapMeasuresImageFilter, which is what we do in the utility function compute_overlap_measures.
Computing the boundary distance scores is slightly more involved. This is a two step process, where first an unsigned distance map from the reference data is generated. Then for each segmentation, the voxels on its boundary are labeled, LabelContour, and the intensity statistics is computed using the label as a mask and the distance map as the image for the LabelIntensityStatisticsImageFilter. This is implemented in the utility function compute_surface_distance_measures.

R> surface_distance_measures <-as.data.frame(surface_distance_measures) R> surface_distance_measures$rater <-rownames(surface_distance_measures)
It is straightforward to visually compare raters using their overlap scores, as illustrated in the code below, with the results shown in Figure 5.

Cell segmentation
Segmentation of cells in fluorescent microscopy is a relatively common image characterization task with variations that are dependent on the specifics of fluorescent markers for a given experiment. A typical procedure might include: • Histogram-based threshold estimation to produce a binary image.
• Cell splitting (separating touching cells) using distance transforms and watershed transform.
• Refinement of initial segmentation using information from other channels.
This example demonstrates the procedure on a 3 channel fluorescent microscopy image. The blue channel is a DNA marker (DAPI) that indicates all cells, the red channel is a marker of cell death (Ph3) while the green channel is a marker of cell proliferation (Ki67). A typical experiment might count the number of cells and measure size in the different states, where states are determined by presence of Ph3 and Ki67, various times after treatment with a drug candidate.

Cell segmentation and splitting
Histogram-based threshold estimation is performed by the segChannel function, listed below. It applies light smoothing followed by the Li threshold estimator (Li and Tam 1998), one of a range of threshold estimation options available in SimpleITK. A cell splitting procedure using distance transforms and a marker-based watershed (implemented by segBlobs, also listed below) was then applied to the resulting mask. Distance transforms replace each foreground pixel with the distance to the closest background pixel, producing a cone-shaped brightness profile for each circular object. Touching cells can then be separated using the peaks of the cones as markers in a watershed transform. A marker image is created by identifying peaks in the distance transform and applying a connected-component labeling. The inverted distance transform is used as the control image for the watershed transform.
The original image, provided in Nowell (2015), is illustrated in Figure 6 and processing stages are illustrated for an image subset (

Refinement of segmentation
The Ph3 and Ki67 stains do not mark the entire cell nuclei. In some cases the marker is much smaller than the nucleus and there may be multiple markers per nucleus, leading to errors in counts or areas. We can, however, use the segmentation results from these channels to perform a geodesic reconstruction based on the DAPI segmentation results to reliably segment cells with those markers. The steps below mask the Ph3 and Ki67 segmentations using the DAPI segmentation and then apply a geodesic reconstruction that "grows" the Ph3/Ki67 cell segmentation to the size of the DAPI marked cell. Figure 8 illustrates the initial segmentation of the Ph3 and Ki67 channels while Figure 9 illustrates the results for the Ph3 channel.   R> dapi.mask <-dapi.cells$labels != 0 R> ph3.markers <-ph3.cells$thresh * dapi.mask R> Ki67.markers <-Ki67.cells$thresh * dapi.mask    Using this data we can begin to visualize the properties of the cell population -for example distributions of cell areas - Figure 10. Note that the cell measures include information about which cells touch the image boundary, allowing easy removal to avoid biasing area statistics.

Discussion and conclusions
There are a large number of R packages with imaging capabilities, many of which are discussed in the CRAN task view on medical imaging (Whitcher 2018) and the recent special issue on "Magnetic Resonance Imaging in R" of this journal (Tabelow and Whitcher 2011). Some packages implement IO of standard or specialized image formats: tiff ( Other packages focus on specific techniques, such as texture analysis (radiomics; Carlson 2016), smoothing (adimpro; Tabelow and Polzehl 2016a), boundary detection (BayesBD; Li 2017), and registration (NiftiReg; Modat et al. 2010;and antsR;Avants, Kandel, Duda, Cook, and Tustison 2015).
Finally, some packages provide interfaces to external imaging applications, such as fslr (Muschelli 2018) and RImageJ (Francois, Grosjean, and Murrell 2015), interfacing to FSL and ImageJ (Schneider, Rasband, and Eliceiri 2012;Schindelin et al. 2012) respectively. R is acting as an advanced shell to run executables from those packages.
Although many of the packages above may have functionality that is applicable to multiple domains, they are not designed to be general purpose packages. Two exceptions are imager (Barthelme 2017), which utilizes the CImg library to provide a range of filtering, thresholding and warping functions and EBImage (Pau, Fuchs, Sklyar, Boutros, and Huber 2010) which provides similar capabilities for cell imaging in microscopy. Both use R arrays to represent images. SimpleITK offers a much greater range of computational functions than either of these packages.
The SimpleITK package introduced in this article supports IO for a wide variety of image formats and provides a broad range of computational components to perform filtering, segmentation and registration. It uses the image data structures from ITK, which include a complete set of meta-data describing image and voxel geometry in world coordinates, and processing classes, which have been widely used for many years and have solid software engineering support to provide long-term maintainability.
SimpleITK provides an open source solution to complex image analysis tasks with facilities that rival many commercial offerings. The option of using SimpleITK with interpreter-based environments such as R allows a researcher or developer to quickly explore many combinations of computational strategies when working with images.
To follow the ongoing toolkit development go to: https://github.com/SimpleITK/ SimpleITK. We hope that SimpleITK will be useful to anyone faced with the task of obtaining insights from images.