Images and pixels with SimpleITK from R¶
The R notebooks follow the general structure of the python equivalents.
A few comments to get started:
- R will display a variable when you type the name. If the variable is a SimpleITK image then, by default, R/SimpleITK will attempt to display it with ImageJ. This isn't much good for notebooks or knitted documents, so the "show" method has been replaced by a function using internal R graphics. This isn't necessary for interactive sessions.
- The documentation doesn't exist in the R package yet - be prepared to do some detective work.
- C++ enumerated types are wrapped using strings, rather than integers.
- Error messages are can be very cryptic - it helps to get familiar with the R way of querying methods and objects.
- The object oriented style used by swig is one of the less common ones used in R packages.
- C++ vectors are automatically converted to R vectors.
Let's start by creating an image and accessing information about it.
# Load the SimpleITK library
library(SimpleITK)
imA = Image(256, 128, 64, "sitkInt16")
imB = SimpleITK::Image(64, 64, "sitkFloat32")
imC = Image(c(32, 32), "sitkUInt32")
imRGB = Image(c(128,128), "sitkVectorUInt8", 3)
We've created 3 different (empty) images using slightly differing syntax to specify the dimensions. We've also specified the pixel type. The range of available pixel types is:
getAnywhere(".__E___itk__simple__PixelIDValueEnum")
It isn't usually necessary to access the enumeration classes directly, but it can be handy to know how when debugging.
Lets start querying the images we've created. The $ operator is used to access object methods:
imA$GetSize()
imC$GetSpacing()
We can change the spacing (voxel size):
imC$SetSpacing(c(2, 0.25))
imC$GetSpacing()
There are shortcuts for accessing a single dimension:
imC$GetWidth()
imA$GetHeight()
imA$GetDepth()
imC$GetDepth()
The depth refers to the z dimension. A colour or vector image has an additional "dimension":
imRGB$GetNumberOfComponentsPerPixel()
The list of available methods can be retrieved using the following. There are many functions for accessing pixels in a type dependent way, which don't need to be accessed directly.
getMethod('$', class(imA))
Array notation for images¶
R, like python and MATLAB, has a flexible notation for accessing elements of arrays. SimpleITK adopts the same notation for accessing pixels (or sub images, which will be covered in a subsequent notebook). First, lets get set up so that we can display an image in the notebook.
# override show function
# The example below will respect image voxel sizes
## This is a hack to allow global setting of the displayed image width
Dwidth <- grid::unit(10, "cm")
setMethod('show', '_p_itk__simple__Image',
function(object)
{
a <- t(as.array(object))
rg <- range(a)
A <- (a-rg[1])/(rg[2]-rg[1])
dd <- dim(a)
sp <- object$GetSpacing()
sz <- object$GetSize()
worlddim <- sp*sz
worlddim <- worlddim/worlddim[1]
W <- Dwidth
H <- Dwidth * worlddim[2]
grid::grid.raster(A,default.units="mm", width=W, height=H)
}
)
options(repr.plot.height = 5, repr.plot.width=5)
Typing the name of an image variable will now result in the image being displayed using an R graphics device. Lets load up an image from the package and test:
Now we source the code to allow us to fetch images used in examples. This code will download the images if they aren't already on the system.
source("downloaddata.R")
cthead <- ReadImage(fetch_data("cthead1.png"))
cthead
Now we will check the value of a single pixel:
cthead[146, 45]
Pixel access order and conversion to arrays¶
The rule for index ordering in R/SimpleITK is most rapidly changing index first. For images this is x (horizontal), then y (vertical), then z. For R arrays this is row, then column. The conversion between images and arrays maintains the ordering of speed of index change, so the most rapidly changing index is always the first element:
imA$GetSize()
imA.arr <- as.array(imA)
dim(imA.arr)
The array and image dimensions are the same, however if they are displayed they are likely to appear different as R interprets the first index as rows. Also note that traditional graphing tools are likely to interpret vertical direction differently: For example, lets plot the cthead image - the result is flipped vertically, because the graph origin is bottom left, and rows and columns are swapped. The show method we use in this notebook has a transpose operation to reverse the row/column swap.
image(as.array(cthead))
