ANTsR.Rmd
“A small leak will sink a great ship.” (folk wisdom)
The ANTsR package interfaces state of the art image
processing with R statistical methods. The project grew out of
the need, at University of Pennsylvania, to develop large-scale
analytics pipelines that track provenance from scanner to scientific
study. ANTsR achieves this by wrapping an ANTs and ITK C++ core
via Rcpp
(Eddelbuettel
2013).
ITK is a templated C++ framework with I/O and support for arbitrary image types (usually 2, 3 or 4 dimensions) as well as surface representations. ANTs, built on ITK, focuses on multivariate image matching and segmentation as well as geometric (even high-dimensional) image transformation. Both tools are deeply validated and widely used.
Together, these tools allow powerful image manipulations. However, they lack a true statistical back-end. Historically, statistical software was not amenable to direct manipulation of multiple dimensional images. This led to “in-house” statistical programming or, perhaps worse (for science), reliance on closed source commercial software. Given the increasing popularity of R and prominence of quantitative imaging, it is natural that R should have a package focused on biological or medical image analysis.
This package integrates several frameworks for extracting quantitative information from images and mapping images into reference coordinate systems. Human brain mapping studies have long relied on Talairach-Tournoux and related coordinate systems (TALAIRACH and TOURNOUX 1958). Similar standardized localization is becoming more common within non-human studies (Johnson et al. 2010; Majka et al. 2013). Atlases of other organ systems are also emerging and being applied clinically (de Marvao et al. 2014). This class of methods relies on image transformation and image segmentation as an aid to the ultimate goal of quantifying variability within and across populations. Longer term, such methods will be critical to individualized patient care and other translational applications.
Here, we provide an overview of the methods available within ANTsR.
core image processing and I/O: ITK (B. B. Avants, Tustison, et al. 2014);
registration and utilities for image processing: ANTs mappings (Tustison, Cook, et al. 2014) and feature extraction (Tustison, Shrinidhi, et al. 2014);
dimensionality reduction: Eigenanatomy (Dhillon et al. 2014) and SCCAN (B. B. Avants, Libon, et al. 2014);
methods for ASL-based cerebral blood flow quantification (Kandel et al. 2015);
neighborhood representations of images that enable rich statistical models (Kandel et al. 2015)
core statistics and temporal filtering via R packages that is amenable to BOLD image processing
In combination, these tools enable one to go from near-raw medical imaging data to a fully reproducible scientific publication (B. B. et al. Avants 2015).
This package uses an antsImage
S4 class to hold pointers
to ITK images. We convert antsImage
objects to R
objects before passing them to R statistical methods. E.g. we
convert a scalar image to a vector, a collection of scalar
images to a matrix or a time series image to a matrix. Currently,
ANTsR does not explicitly represent images with vector valued
voxels (e.g. tensor or warp images) although these may be supported in
the future in a way that is similar to our current support for time
series images. The large majority of images employed within
ANTsR are of 2, 3 or 4 dimensions with float
pixel
types. This information is stored within the antsImage
class.
A few example images are built into ANTsR, but more can be
downloaded. See ?getANTsRData
.
img<-antsImageRead( getANTsRData("r16"), 2 ) # built in image
img
## antsImage
## Pixel Type : float
## Components Per Pixel: 1
## Dimensions : 256x256
## Voxel Spacing : 1x1
## Origin : 0 0
## Direction : 1 0 0 1
## Filename : /Users/stnava/Library/R/arm64/4.3/library/ANTsRCore/extdata/r16slice.jpg
Take a quick look at the image.
ANTsR includes:
An organizational system such that relatively small scripts may implement full studies
Implementation of foundational methods
compcor
and
*DenoiseR
Reference simulation data and examples distributed with the package
Interpretation of results
Openness and reproducibility
In total, ANTsR is a rigorous framework upon which one may
build customized statistical implementations appropriate for large-scale
functional, structural or combined functional and structural image
analyses. Because much of the core is implemented with C++, the
framework also remains efficient. Finally, note that
Rscript
allows one to send ANTsR scripts to
clusters and take advantage of distributed computing resources.
Here, we quickly summarize ANTsR functionality and useful tools.
The travis build system
We test ANTsR regularly. The status of the build (and an expected build result) can be seen here: . Take a look at the detailed log to see what one might expect if building ANTsR from source.
Image input and output
If nothing else, ANTsR makes it easy to read and write
(medical) images and to map them into a format compatible with
R. Formats we frequently use include jpg, tiff, mha, nii.gz and
nrrd. However, only the last three have a proper physical space
representation necessary for mapping. Below is an example of how we
access this type of image and see its geometry. Check
antsImageWrite
for the primary supported I/O.
mnifilename<-getANTsRData("r27")
img<-antsImageRead(mnifilename, pixeltype="unsigned char")
img
## antsImage
## Pixel Type : unsigned char
## Components Per Pixel: 1
## Dimensions : 256x256
## Voxel Spacing : 1x1
## Origin : 0 0
## Direction : 1 0 0 1
## Filename : /Users/stnava/Library/R/arm64/4.3/library/ANTsRCore/extdata/r27slice.jpg
# retval<-antsImageWrite(img,mnifilename)
antsGetSpacing(img)
## [1] 1 1
antsGetDirection(img)
## [,1] [,2]
## [1,] 1 0
## [2,] 0 1
antsGetOrigin(img)
## [1] 0 0
print(img[120,122]) # same type of thing in 3 or 4D
## [,1]
## [1,] 182
## [1] 255
Index an image with a label
Often, you would like to summarize or extract information from within a known region of an image with arbitrary shape but within a given intensity “zone”. We simulate this below and show a few accessors and type conversions.
gaussimg<-array( data=rnorm(125), dim=c(5,5,5) )
arrayimg<-array( data=(1:125), dim=c(5,5,5) )
img<-as.antsImage( arrayimg )
print( max(img) )
## [1] 125
## [1] 88
# if using SUBSET using an antsImage, you must be explicit
sub = as.array(img >= 50) > 0
print( mean( gaussimg[ sub ]) )
## [1] 0.1308089
Convert a 4D image to a matrix
Four dimensional images are generated and used in the same way. One can easily transform from 4D image to matrix and back.
## [1] 5 5 5 10
avg3d<-ANTsR::getAverageOfTimeSeries( gaussimg )
mask <- avg3d < 0.25
gmat<-timeseries2matrix( gaussimg, mask )
print(dim(gmat))
## [1] 10 103
If one has a mask, then one can use makeImage
to
generate a new image from a scalar or vector.
newimg<-makeImage( mask, mean(avg3d) ) # from scalar
newimg<-makeImage( mask, colMeans(gmat) ) # from vector
Convert a list of images to a matrix
Often, one has several scalar images that need to be accumulated for statistical processing. Here, we generate a simulated set of these images and then proceed to smooth them, store them in a list and convert them to a matrix after extracting the information of each image within a data-driven mask.
nimages<-100
ilist<-list()
for ( i in 1:nimages )
{
simimg<-makeImage( c(50,50) , rnorm(2500) )
simimg<-smoothImage(simimg,1.5)
ilist[[ i ]] = simimg
}
# get a mask from the first image
mask<-getMask( ilist[[1]],
lowThresh=mean(ilist[[1]]), cleanup=TRUE )
mat<-imageListToMatrix( ilist, mask )
## | | | 0% | |= | 1% | |= | 2% | |== | 3% | |=== | 4% | |==== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 16% | |============ | 17% | |============= | 18% | |============= | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 23% | |================= | 24% | |================== | 25% | |================== | 26% | |=================== | 27% | |==================== | 28% | |==================== | 29% | |===================== | 30% | |====================== | 31% | |====================== | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 36% | |========================== | 37% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================= | 41% | |============================= | 42% | |============================== | 43% | |=============================== | 44% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 56% | |======================================== | 57% | |========================================= | 58% | |========================================= | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================= | 64% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 83% | |=========================================================== | 84% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%
## [1] 100 362
Once we have a matrix representation of our population, we might run a quick voxel-wise regression within the mask. Then we look at some summary statistics.
mat<-imageListToMatrix( ilist, mask )
## | | | 0% | |= | 1% | |= | 2% | |== | 3% | |=== | 4% | |==== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 16% | |============ | 17% | |============= | 18% | |============= | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 23% | |================= | 24% | |================== | 25% | |================== | 26% | |=================== | 27% | |==================== | 28% | |==================== | 29% | |===================== | 30% | |====================== | 31% | |====================== | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 36% | |========================== | 37% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================= | 41% | |============================= | 42% | |============================== | 43% | |=============================== | 44% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 56% | |======================================== | 57% | |========================================= | 58% | |========================================= | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================= | 64% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 83% | |=========================================================== | 84% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%
age<-rnorm( nrow(mat) ) # simulated age
gender<-rep( c("F","M"), nrow(mat)/2 ) # simulated gender
# this creates "real" but noisy effects to detect
mat<-mat*(age^2+rnorm(nrow(mat)))
mydf = data.frame( age = age, gender=factor(gender ))
mdl<-lm( mat ~ age + gender, data=mydf )
mdli<-bigLMStats( mdl, 1.e-4 )
print(names(mdli))
## [1] "fstat" "pval.model" "beta" "beta.std" "beta.t"
## [6] "beta.pval"
## [1] "age" "genderM"
## [1] "age 0.0570923169814"
## [1] "gen 1"
Write out a statistical map
We might also write out the images so that we can save them for later or look at them with other software.
agebetas<-makeImage( mask , mdli$beta.t[1,] )
returnval<-antsImageWrite( agebetas, tempfile(fileext ='.nii.gz') )
We achieve quantification in biological or medical imaging by using prior knowledge about the image content.
Segmentation
In segmentation, we assume the image has a known set of tissues, organs etc. Here, we assume 3 tissues exist and use a classic k-means model with MRF penalty (B. B. Avants et al. 2011). Note that we also bias correct the image to help it match our model (Tustison et al. 2010).
fi<-antsImageRead( getANTsRData("r16") ,2)
fi<-n3BiasFieldCorrection(fi,2)
seg<-kmeansSegmentation( fi, 3 )
invisible(plot(seg$segmentation))
If you like segmentation, also look at rfSegmentation
and atropos
.
Registration
In registration, we assume the image can be mapped to some canonical shape or example, i.e. an atlas. Or to another individual. ANTsR provides a simple wrapper for SyN image registration (Tustison and Avants 2013),
mi<-antsImageRead( getANTsRData("r64") ,2)
mytx<-antsRegistration(fixed=fi , moving=mi ,
typeofTransform = c('SyN'))
regresult<-iMath(mytx$warpedmovout,"Normalize")
fiedge<-iMath(fi,"Canny",1,5,12)
invisible(plot(regresult, list(fiedge), window.overlay=c(0.5,1)) )
while invariantImageSimilarity
provides powerful
multi-start search for lower dimensional affine registrations.
Deformable image registration results in a voxel-wise map of the contraction and expansion of the moving image (after affine transformation) that is needed to map to the fixed image. This deformation gradient is colloquially known as “the jacobian”.
jac<-createJacobianDeterminantImage(fi,mytx$fwdtransforms[[1]],1)
invisible(plot(jac))
Above, we compute and plot the image of the log-jacobian. This mapping is a useful summary measurement for morphometry Kim et al. (2008).
Registration and segmentation
Registration and segmentation are often applied jointly or
iteratively to maximize some criterion. See the example in
jointIntensityFusion
for one such case (Wang and Yushkevich 2013).
Neighborhood operations
Basic I/O and management of images as vectors is critical. However, there is additional information that can be gained by representing an image and its neighborhood information. ANTsR represents image neighborhoods, which capture shape and texture, as a matrix. Here, extract an image neighborhood matrix representation such that we may analyze it at a given scale.
mnit<-getANTsRData("r16")
mnit<-antsImageRead( mnit )
mnit <- resampleImage( mnit , rep(4, mnit@dimension) ) # downsample
mask2<-getMask(mnit,lowThresh=mean(mnit),cleanup=TRUE)
radius <- rep(2,mnit@dimension)
mat2<-getNeighborhoodInMask(mnit, mask2, radius,
physical.coordinates = FALSE,
boundary.condition = "mean" )
print(dim(mat2))
## [1] 25 1113
The variable mat2
has size determined by the
neighborhood radius (here, 5) and the number of non-zero voxels in the
mask. The boundary.condition
says how to treat data that is
outside of the mask or the image boundaries. This example replaces
missing data with the mean in-mask value of the local neighborhood.
Other useful tools in ANTsR include iMath
,
thresholdImage
, quantifyCBF
,
preprocessfMRI
, aslPerfusion
,
computeDVARS
, getROIValues
,
hemodynamicRF
, makeGraph
,
matrixToImages
, rfSegmentation
,
antsRegistration
, plotPrettyGraph
,
plotBasicNetwork
, getTemplateCoordinates
,
antsSet*
.
Several image mathematics operations (like ImageMath
in
ANTs) are accessible too via iMath
.
ANTsR also provides AAL label (Tzourio-Mazoyer et al. 2002) names via:
data(aal,package='ANTsR')
labs<-1:90
with cortical labs defined by labs
. The DKT atlas labels
are similarly summarized in DesikanKillianyTourville
(Klein and Tourville 2012).
An example BOLD correlation matrix is available in
bold_correlation_matrix
. This can be used to try out
makeGraph
and related functions.
The basic plot
function is implemented for the
antsImage
class. It can show 2 or 3D data with color
overlays, the latter of which can display multiple slices side by side.
Several color choices are available for the overlays.
For 3D images, see renderSurfaceFunction
and
plotBasicNetwork
for rgl
and
misc3d
based interactive surface and network plots. Another
such example is in visualizeBlob
. These are too
long-running to compile into the vignette but the help examples for
these functions will allow you to see their results.
A good visualization alternative outside of ANTsR is antsSurf.
Good approaches exist in ANTsR for preprocessing BOLD data.
These yield both motion matrices and relevant summary measurements such
as FD and DVARS. See ?preprocessfMRI
for a simplified
utility function. This function could be used on each run of an
experiment and the results stored in organized fashion for later
use.
Motion correction
To motion correct your data, one might run:
# get an average image
averageImage <- getAverageOfTimeSeries( boldImage )
motionCorrectionResults <- antsMotionCalculation( boldImage,
fixed = averageImage )
A moreaccurate
flag should be set to 1
or
2
for usable (not test) results. FD and DVARS are returned
which may be used to summarize motion. One might also get this data from
preprocessfMRI
which also provides denoising options based
on data-driven methods including frequency filtering.
For more fMRI focused tools, see RKRNS and its github site github RKRNS.
Images often have many voxels (-voxels) and, in medical applications, this means that or even , where is the number of subjects. Therefore, we often want to “intelligently” reduce the dimensionality of the data. We favor methods related to PCA and CCA but have a few ICA related tools too.
Eigenanatomy & SCCAN
Our sparse and geometrically constrained dimensionality reduction methods seek to both explain variance and also yield interpretable, spatially localized pseudo-eigenvectors Cook et al. (2014). This is the point of “eigenanatomy” which is a variation of sparse PCA that uses (optionally) biologically-motivated smoothness, locality or sparsity constraints.
# assume you ran the population example above
eanat<-sparseDecom( mat, mask, 0.2, 5, cthresh=2, its=2 )
eanatimages = matrixToImages( eanat$eig, mask )
eseg<-eigSeg(mask, eanatimages ,F)
jeanat<-joinEigenanatomy(mat, mask, eanatimages, c(0.1))
eseg2<-eigSeg(mask,jeanat$fusedlist,F)
The parameters for the example above are set for fast processing. You can see our paper for some theory on these methods (Kandel et al. 2014). A more realistic study setup would be
eanat<-sparseDecom( inmatrix=mat, inmask=famask, nvecs=50,
sparseness=0.005, cthresh=500, its=5, mycoption=0 )
jeanat<-joinEigenanatomy( mat , famask, eanat$eig,
c(1:20)/100.0 , joinMethod='multilevel' )
useeig<-eanat$eig
useeig<-jeanat$fusedlist
avgmat<-abs(imageListToMatrix( useeig , famask ))
avgmat<-avgmat/rowSums(abs(avgmat))
imgmat<-( mat %*% t(avgmat) )
The imgmat
variable would be your summary predictors
entered into lm
or randomForest
.
More information is available within the examples that can be seen
within the help for sparseDecom
, sparseDecom2
and the helper function initializeEigenanatomy
.
Sparse canonical correlation analysis
CCA maximizes
where
are as above and
and
are similarly defined. CCA optimizes the matrices
operating on
to find a low-dimensional representation of the data pair
in which correlation is maximal. Following ideas outlined in Dhillon et al. (2014) and B. B. Avants, Libon, et al. (2014), this method
can be extended with sparsity constraints that yield rows of
with a controllable number of non-zero entries. See the sccan tutorial and
sparseDecom2
for more information.
With the current ANTsR, one may:
Exploit ANTs and ITK functionality within R
Leverage R functionality to help understand and interpret imaging data
Use feature selection based on various filtering strategies in
iMath
and elsewhere (e.g
segmentShapeFromImage
)
Employ dimensionality reduction through eigenanatomy or SCCAN with a variety of incarnations, some of which are similar to ICA
Use relatively few interpretable and low-dimensional predictors derived from high-dimensional data.
Interpret multivariate results intuitively when used in combination with standard R visualization.
See ANTsR for all source code and documentation and RKRNS-talk for html slides that discuss extensions to BOLD decoding.
Enjoy and please refer issues to ANTsR issues.