-
Notifications
You must be signed in to change notification settings - Fork 35
MNI Coordinates in ANTsR (and ANTs)
Let us suppose we want the MNI coordinate of posterior cingulate cortex. We choose something from the literature which gives 0.0, 53.0, 26.0
. In ANTsR
, these coordinates work perfectly well. Choose another example and the result looks wrong. The difference relates to the underlying coordinate system of the neuroimage representation, e.g. whether the data is in LPS or another coordinate system. The difference is often in the sign of the y
coordinate but we recommend checking explicitly in your own data. Here is a way to check the anatomical location of coordinates in ANTsR
.
# create a sphere image
if ( !exists("mask") )
mask = antsImageRead( getANTsRData( "mni" ) ) %>% resampleImage( c( 4, 4, 4 ) )
pccMask = mask * 0
pccCoord = as.numeric( c( 0.0, 53.0, 26.0 ) ) # there is variability in reported MNI coords of PCC
idx = antsTransformPhysicalPointToIndex( mask, pccCoord )
rad = 6.0 # millimeters
n = ceiling( rad / antsGetSpacing( mask ) )
for (i in c(-n[1]:n[1]) ) {
for (j in c(-n[2]:n[2])) {
for (k in c(-n[3]:n[3])) {
local = idx + c(i,j,k)
localpt = as.numeric( antsTransformIndexToPhysicalPoint( mask, local ) )
dist = sqrt( sum( (localpt - pccCoord)*(localpt - pccCoord) ))
inImage = (prod(idx <= dim(mask)) == 1) && (length(which(idx < 1)) == 0 )
if ( (dist <= rad) && (inImage == TRUE ) ) {
pccMask[ local[1], local[2], local[3] ] = 1
}
}
}
}
plot( mask, pccMask, axis = 3, nslices = 30, ncolumns = 10 )
Now here is the brief version using an existing ANTsR function.
# create a sphere image
if ( !exists("mask") )
mask = antsImageRead( getANTsRData( "mni" ) ) %>% resampleImage( c( 4, 4, 4 ) )
pccMask = makePointsImage( matrix( c(0, 52, 34), nrow=1 ), mask, 10 )
plot( mask, pccMask, axis = 3, nslices = 30, ncolumns = 10 )
Now let us suppose we want to map MNI points to our own template. We use an invertible map between 2 domains and apply to a point basis.
fixed <- antsImageRead( getANTsRData("ch2") )
moving <- antsImageRead( getANTsRData("mni") )
mytx <- antsRegistration( fixed=fixed, moving=moving, typeofTransform = c("SyN") )
data( "powers_areal_mni_itk", package = "ANTsR", envir = environment() )
coords = powers_areal_mni_itk[,1:3]
ch2reg = antsRegistration( fixed, moving, typeofTransform = "SyN" )
coordsw <- antsApplyTransformsToPoints( dim=3, points=coords,
transformlist=ch2reg$fwdtransforms, whichtoinvert=c(FALSE,FALSE) )
ptrd = 3
powersLabels = makePointsImage( coordsw, moving, radius = pard )
plot( moving, powersLabels, axis=3, nslices=30 )