#' Cortical and deep gray matter labeling using Desikan-Killiany-Tourville
#'
#' Perform DKT labeling using deep learning
#'
#' The labeling is as follows:
#' \itemize{
#' \item{Label 0:}{background}
#'
#' Inner labels:
#' \item{Label 4:}{left lateral ventricle}
#' \item{Label 5:}{left inferior lateral ventricle}
#' \item{Label 6:}{left cerebellem exterior}
#' \item{Label 7:}{left cerebellum white matter}
#' \item{Label 10:}{left thalamus proper}
#' \item{Label 11:}{left caudate}
#' \item{Label 12:}{left putamen}
#' \item{Label 13:}{left pallidium}
#' \item{Label 14:}{3rd ventricle}
#' \item{Label 15:}{4th ventricle}
#' \item{Label 16:}{brain stem}
#' \item{Label 17:}{left hippocampus}
#' \item{Label 18:}{left amygdala}
#' \item{Label 24:}{CSF}
#' \item{Label 25:}{left lesion}
#' \item{Label 26:}{left accumbens area}
#' \item{Label 28:}{left ventral DC}
#' \item{Label 30:}{left vessel}
#' \item{Label 43:}{right lateral ventricle}
#' \item{Label 44:}{right inferior lateral ventricle}
#' \item{Label 45:}{right cerebellum exterior}
#' \item{Label 46:}{right cerebellum white matter}
#' \item{Label 49:}{right thalamus proper}
#' \item{Label 50:}{right caudate}
#' \item{Label 51:}{right putamen}
#' \item{Label 52:}{right palladium}
#' \item{Label 53:}{right hippocampus}
#' \item{Label 54:}{right amygdala}
#' \item{Label 57:}{right lesion}
#' \item{Label 58:}{right accumbens area}
#' \item{Label 60:}{right ventral DC}
#' \item{Label 62:}{right vessel}
#' \item{Label 72:}{5th ventricle}
#' \item{Label 85:}{optic chasm}
#' \item{Label 91:}{left basal forebrain}
#' \item{Label 92:}{right basal forebrain}
#' \item{Label 630:}{cerebellar vermal lobules I-V}
#' \item{Label 631:}{cerebellar vermal lobules VI-VII}
#' \item{Label 632:}{cerebellar vermal lobules VIII-X}
#'
#' Outer labels:
#' \item{Label 1002:}{left caudal anterior cingulate}
#' \item{Label 1003:}{left caudal middle frontal}
#' \item{Label 1005:}{left cuneus}
#' \item{Label 1006:}{left entorhinal}
#' \item{Label 1007:}{left fusiform}
#' \item{Label 1008:}{left inferior parietal}
#' \item{Label 1009:}{left inferior temporal}
#' \item{Label 1010:}{left isthmus cingulate}
#' \item{Label 1011:}{left lateral occipital}
#' \item{Label 1012:}{left lateral orbitofrontal}
#' \item{Label 1013:}{left lingual}
#' \item{Label 1014:}{left medial orbitofrontal}
#' \item{Label 1015:}{left middle temporal}
#' \item{Label 1016:}{left parahippocampal}
#' \item{Label 1017:}{left paracentral}
#' \item{Label 1018:}{left pars opercularis}
#' \item{Label 1019:}{left pars orbitalis}
#' \item{Label 1020:}{left pars triangularis}
#' \item{Label 1021:}{left pericalcarine}
#' \item{Label 1022:}{left postcentral}
#' \item{Label 1023:}{left posterior cingulate}
#' \item{Label 1024:}{left precentral}
#' \item{Label 1025:}{left precuneus}
#' \item{Label 1026:}{left rostral anterior cingulate}
#' \item{Label 1027:}{left rostral middle frontal}
#' \item{Label 1028:}{left superior frontal}
#' \item{Label 1029:}{left superior parietal}
#' \item{Label 1030:}{left superior temporal}
#' \item{Label 1031:}{left supramarginal}
#' \item{Label 1034:}{left transverse temporal}
#' \item{Label 1035:}{left insula}
#' \item{Label 2002:}{right caudal anterior cingulate}
#' \item{Label 2003:}{right caudal middle frontal}
#' \item{Label 2005:}{right cuneus}
#' \item{Label 2006:}{right entorhinal}
#' \item{Label 2007:}{right fusiform}
#' \item{Label 2008:}{right inferior parietal}
#' \item{Label 2009:}{right inferior temporal}
#' \item{Label 2010:}{right isthmus cingulate}
#' \item{Label 2011:}{right lateral occipital}
#' \item{Label 2012:}{right lateral orbitofrontal}
#' \item{Label 2013:}{right lingual}
#' \item{Label 2014:}{right medial orbitofrontal}
#' \item{Label 2015:}{right middle temporal}
#' \item{Label 2016:}{right parahippocampal}
#' \item{Label 2017:}{right paracentral}
#' \item{Label 2018:}{right pars opercularis}
#' \item{Label 2019:}{right pars orbitalis}
#' \item{Label 2020:}{right pars triangularis}
#' \item{Label 2021:}{right pericalcarine}
#' \item{Label 2022:}{right postcentral}
#' \item{Label 2023:}{right posterior cingulate}
#' \item{Label 2024:}{right precentral}
#' \item{Label 2025:}{right precuneus}
#' \item{Label 2026:}{right rostral anterior cingulate}
#' \item{Label 2027:}{right rostral middle frontal}
#' \item{Label 2028:}{right superior frontal}
#' \item{Label 2029:}{right superior parietal}
#' \item{Label 2030:}{right superior temporal}
#' \item{Label 2031:}{right supramarginal}
#' \item{Label 2034:}{right transverse temporal}
#' \item{Label 2035:}{right insula}
#' }
#'
#' Performing the lobar parcellation is based on the FreeSurfer division
#' described here:
#'
#' \url{https://surfer.nmr.mgh.harvard.edu/fswiki/CorticalParcellation}
#'
#' \itemize{
#' Frontal lobe:
#' \item{Label 1002:}{left caudal anterior cingulate}
#' \item{Label 1003:}{left caudal middle frontal}
#' \item{Label 1012:}{left lateral orbitofrontal}
#' \item{Label 1014:}{left medial orbitofrontal}
#' \item{Label 1017:}{left paracentral}
#' \item{Label 1018:}{left pars opercularis}
#' \item{Label 1019:}{left pars orbitalis}
#' \item{Label 1020:}{left pars triangularis}
#' \item{Label 1024:}{left precentral}
#' \item{Label 1026:}{left rostral anterior cingulate}
#' \item{Label 1027:}{left rostral middle frontal}
#' \item{Label 1028:}{left superior frontal}
#' \item{Label 2002:}{right caudal anterior cingulate}
#' \item{Label 2003:}{right caudal middle frontal}
#' \item{Label 2012:}{right lateral orbitofrontal}
#' \item{Label 2014:}{right medial orbitofrontal}
#' \item{Label 2017:}{right paracentral}
#' \item{Label 2018:}{right pars opercularis}
#' \item{Label 2019:}{right pars orbitalis}
#' \item{Label 2020:}{right pars triangularis}
#' \item{Label 2024:}{right precentral}
#' \item{Label 2026:}{right rostral anterior cingulate}
#' \item{Label 2027:}{right rostral middle frontal}
#' \item{Label 2028:}{right superior frontal}
#'
#' Parietal:
#' \item{Label 1008:}{left inferior parietal}
#' \item{Label 1010:}{left isthmus cingulate}
#' \item{Label 1022:}{left postcentral}
#' \item{Label 1023:}{left posterior cingulate}
#' \item{Label 1025:}{left precuneus}
#' \item{Label 1029:}{left superior parietal}
#' \item{Label 1031:}{left supramarginal}
#' \item{Label 2008:}{right inferior parietal}
#' \item{Label 2010:}{right isthmus cingulate}
#' \item{Label 2022:}{right postcentral}
#' \item{Label 2023:}{right posterior cingulate}
#' \item{Label 2025:}{right precuneus}
#' \item{Label 2029:}{right superior parietal}
#' \item{Label 2031:}{right supramarginal}
#'
#' Temporal:
#' \item{Label 1006:}{left entorhinal}
#' \item{Label 1007:}{left fusiform}
#' \item{Label 1009:}{left inferior temporal}
#' \item{Label 1015:}{left middle temporal}
#' \item{Label 1016:}{left parahippocampal}
#' \item{Label 1030:}{left superior temporal}
#' \item{Label 1034:}{left transverse temporal}
#' \item{Label 2006:}{right entorhinal}
#' \item{Label 2007:}{right fusiform}
#' \item{Label 2009:}{right inferior temporal}
#' \item{Label 2015:}{right middle temporal}
#' \item{Label 2016:}{right parahippocampal}
#' \item{Label 2030:}{right superior temporal}
#' \item{Label 2034:}{right transverse temporal}
#'
#' Occipital:
#' \item{Label 1005:}{left cuneus}
#' \item{Label 1011:}{left lateral occipital}
#' \item{Label 1013:}{left lingual}
#' \item{Label 1021:}{left pericalcarine}
#' \item{Label 2005:}{right cuneus}
#' \item{Label 2011:}{right lateral occipital}
#' \item{Label 2013:}{right lingual}
#' \item{Label 2021:}{right pericalcarine}
#'
#' Other outer labels:
#' \item{Label 1035:}{left insula}
#' \item{Label 2035:}{right insula}
#' }
#'
#' Preprocessing on the training data consisted of:
#' * n4 bias correction,
#' * denoising,
#' * brain extraction, and
#' * affine registration to MNI.
#' The input T1 should undergo the same steps. If the input T1 is the raw
#' T1, these steps can be performed by the internal preprocessing, i.e. set
#' \code{doPreprocessing = TRUE}
#'
#' @param t1 raw or preprocessed 3-D T1-weighted brain image.
#' @param doPreprocessing perform preprocessing. See description above.
#' @param returnProbabilityImages whether to return the two sets of probability images
#' for the inner and outer labels.
#' @param doLobarParcellation perform lobar parcellation (also divided by hemisphere).
#' @param antsxnetCacheDirectory destination directory for storing the downloaded
#' template and model weights. Since these can be resused, if
#' \code{is.null(antsxnetCacheDirectory)}, these data will be downloaded to the
#' subdirectory ~/.keras/ANTsXNet/.
#' @param verbose print progress.
#' @return list consisting of the segmentation image and probability images for
#' each label.
#' @author Tustison NJ
#' @examples
#' \dontrun{
#' library( ANTsRNet )
#' library( keras )
#'
#' image <- antsImageRead( "t1.nii.gz" )
#' results <- desikanKillianyTourvilleLabeling( image )
#' }
#' @export
desikanKillianyTourvilleLabeling <- function( t1, doPreprocessing = TRUE,
returnProbabilityImages = FALSE, doLobarParcellation = FALSE,
antsxnetCacheDirectory = NULL, verbose = FALSE )
{
if( t1@dimension != 3 )
{
stop( "Input image dimension must be 3." )
}
################################
#
# Preprocess image
#
################################
t1Preprocessed <- t1
if( doPreprocessing )
{
t1Preprocessing <- preprocessBrainImage( t1,
truncateIntensity = c( 0.01, 0.99 ),
brainExtractionModality = "t1",
template = "croppedMni152",
templateTransformType = "antsRegistrationSyNQuickRepro[a]",
doBiasCorrection = TRUE,
doDenoising = TRUE,
antsxnetCacheDirectory = antsxnetCacheDirectory,
verbose = verbose )
t1Preprocessed <- t1Preprocessing$preprocessedImage * t1Preprocessing$brainMask
}
################################
#
# Download spatial priors for outer model
#
################################
if( verbose )
{
cat( "DesikanKillianyTourville: retrieving label spatial priors.\n" )
}
priorsFileNamePath <- getANTsXNetData( "priorDktLabels",
antsxnetCacheDirectory = antsxnetCacheDirectory )
priorsImageList <- splitNDImageToList( antsImageRead( priorsFileNamePath ) )
################################
#
# Build outer model and load weights
#
################################
templateSize <- c( 96L, 112L, 96L )
labels <- c( 0, 1002:1003, 1005:1031, 1034:1035, 2002:2003, 2005:2031, 2034:2035 )
channelSize <- 1 + length( priorsImageList )
unetModel <- createUnetModel3D( c( templateSize, channelSize ),
numberOfOutputs = length( labels ), mode = 'classification',
numberOfLayers = 4, numberOfFiltersAtBaseLayer = 16, dropoutRate = 0.0,
convolutionKernelSize = c( 3, 3, 3 ), deconvolutionKernelSize = c( 2, 2, 2 ),
weightDecay = 1e-5, additionalOptions = c( "attentionGating" ) )
if( verbose )
{
cat( "DesikanKillianyTourville: retrieving model weights.\n" )
}
weightsFileNamePath <- getPretrainedNetwork( "dktOuterWithSpatialPriors",
antsxnetCacheDirectory = antsxnetCacheDirectory )
load_model_weights_hdf5( unetModel, filepath = weightsFileNamePath )
################################
#
# Do prediction and normalize to native space
#
################################
if( verbose )
{
cat( "DesikanKillianyTourville: outer model prediction.\n" )
}
downsampledImage <- resampleImage( t1Preprocessed, templateSize, useVoxels = TRUE, interpType = 0 )
imageArray <- as.array( downsampledImage )
imageArray <- ( imageArray - mean( imageArray ) ) / sd( imageArray )
batchX <- array( data = 0, dim = c( 1, templateSize, channelSize ) )
batchX[1,,,,1] <- imageArray
for( i in seq.int( length( priorsImageList ) ) )
{
priorImageArray <- as.array( resampleImage( priorsImageList[[i]], templateSize, useVoxels = TRUE, interpType = 0 ) )
batchX[1,,,,i+1] <- priorImageArray
}
predictedData <- unetModel %>% predict( batchX, verbose = verbose )
probabilityImagesList <- decodeUnet( predictedData, downsampledImage )
outerProbabilityImages <- list()
for( i in seq.int( length( probabilityImagesList[[1]] ) ) )
{
resampledImage <- resampleImage( probabilityImagesList[[1]][[i]], dim( t1Preprocessed ), useVoxels = TRUE, interpType = 0 )
if( doPreprocessing )
{
outerProbabilityImages[[i]] <- antsApplyTransforms( fixed = t1, moving = resampledImage,
transformlist = t1Preprocessing$templateTransforms$invtransforms,
whichtoinvert = c( TRUE ), interpolator = "linear", verbose = verbose )
} else {
outerProbabilityImages[[i]] <- resampledImage
}
}
imageMatrix <- imageListToMatrix( outerProbabilityImages, t1 * 0 + 1 )
segmentationMatrix <- matrix( apply( imageMatrix, 2, which.max ), nrow = 1 )
segmentationImage <- matrixToImages( segmentationMatrix, t1 * 0 + 1 )[[1]]
dktLabelImage <- antsImageClone( segmentationImage )
for( i in seq.int( length( labels ) ) )
{
dktLabelImage[( segmentationImage == i )] <- labels[i]
}
################################
#
# Build inner model and load weights
#
################################
templateSize <- c( 160L, 192L, 160L )
labels <- c( 0, 4, 6, 7, 10:18, 24, 26, 28, 30, 43, 44:46, 49:54, 58, 60, 91:92, 630:632 )
unetModel <- createUnetModel3D( c( templateSize, 1 ),
numberOfOutputs = length( labels ), mode = 'classification',
numberOfLayers = 4, numberOfFiltersAtBaseLayer = 8, dropoutRate = 0.0,
convolutionKernelSize = c( 3, 3, 3 ), deconvolutionKernelSize = c( 2, 2, 2 ),
weightDecay = 1e-5, additionalOptions = c( "attentionGating" ) )
if( verbose )
{
cat( "DesikanKillianyTourville: retrieving inner model weights.\n" )
}
weightsFileName <- getPretrainedNetwork( "dktInner",
antsxnetCacheDirectory = antsxnetCacheDirectory )
load_model_weights_hdf5( unetModel, filepath = weightsFileName )
################################
#
# Do inner model prediction and normalize to native space
#
################################
if( verbose )
{
cat( "Inner model prediction.\n" )
}
croppedImage <- cropIndices( t1Preprocessed, c( 13, 15, 1 ), c( 172, 206, 160 ) )
imageArray <- as.array( croppedImage )
batchX <- array( data = imageArray, dim = c( 1, templateSize, 1 ) )
batchX <- ( batchX - mean( batchX ) ) / sd( batchX )
predictedData <- unetModel %>% predict( batchX, verbose = verbose )
probabilityImagesList <- decodeUnet( predictedData, croppedImage )
innerProbabilityImages <- list()
for( i in seq.int( length( probabilityImagesList[[1]] ) ) )
{
if( i > 1 )
{
decroppedImage <- decropImage( probabilityImagesList[[1]][[i]], t1Preprocessed * 0 )
} else {
decroppedImage <- decropImage( probabilityImagesList[[1]][[i]], t1Preprocessed * 0 + 1 )
}
if( doPreprocessing )
{
innerProbabilityImages[[i]] <- antsApplyTransforms( fixed = t1, moving = decroppedImage,
transformlist = t1Preprocessing$templateTransforms$invtransforms,
whichtoinvert = c( TRUE ), interpolator = "linear", verbose = verbose )
} else {
innerProbabilityImages[[i]] <- decroppedImage
}
}
imageMatrix <- imageListToMatrix( innerProbabilityImages, t1 * 0 + 1 )
segmentationMatrix <- matrix( apply( imageMatrix, 2, which.max ), nrow = 1 )
segmentationImage <- matrixToImages( segmentationMatrix, t1 * 0 + 1 )[[1]]
################################
#
# Incorporate the inner model results into the final label image.
# Note that we purposely prioritize the inner label results.
#
################################
for( i in seq.int( length( labels ) ) )
{
if( labels[i] > 0 )
{
dktLabelImage[( segmentationImage == i )] <- labels[i]
}
}
if( doLobarParcellation )
{
if( verbose )
{
cat( "Doing lobar parcellation.\n" )
}
################################
#
# Lobar/hemisphere parcellation
#
################################
# Consolidate lobar cortical labels
if( verbose )
{
cat( " Consolidating cortical labels.\n" )
}
frontalLabels <- c( 1002, 1003, 1012, 1014, 1017:1020, 1024, 1026:1028,
2002, 2003, 2012, 2014, 2017:2020, 2024, 2026:2028 )
parietalLabels <- c( 1008, 1010, 1022, 1023, 1025, 1029, 1031,
2008, 2010, 2022, 2023, 2025, 2029, 2031 )
temporalLabels <- c( 1006, 1007, 1009, 1015, 1016, 1030, 1034,
2006, 2007, 2009, 2015, 2016, 2030, 2034 )
occipitalLabels <- c( 1005, 1011, 1013, 1021,
2005, 2011, 2013, 2021 )
lobarLabels <- list( frontalLabels, parietalLabels, temporalLabels, occipitalLabels )
dktLobes <- antsImageClone( dktLabelImage )
dktLobes[dktLobes < 1000] <- 0
for( i in seq.int( length( lobarLabels ) ) )
{
for( j in seq.int( length( lobarLabels[[i]] ) ) )
{
dktLobes[dktLobes == lobarLabels[[i]][j]] <- i
}
}
dktLobes[dktLobes > length( lobarLabels )] <- 0
sixTissue <- deepAtropos( t1Preprocessed, doPreprocessing = FALSE,
antsxnetCacheDirectory = antsxnetCacheDirectory, verbose = verbose )
atroposSeg <- sixTissue$segmentationImage
if( doPreprocessing )
{
atroposSeg <- antsApplyTransforms( fixed = t1, moving = atroposSeg,
transformlist = t1Preprocessing$templateTransforms$invtransforms,
whichtoinvert = c( TRUE ), interpolator = "genericLabel", verbose = verbose )
}
brainMask <- antsImageClone( atroposSeg )
brainMask[brainMask == 1 | brainMask == 5 | brainMask == 6] <- 0
brainMask <- thresholdImage( brainMask, 0, 0, 0, 1 )
lobarParcellation <- iMath( brainMask, "PropagateLabelsThroughMask", brainMask * dktLobes )
lobarParcellation[atroposSeg == 5] <- 5
lobarParcellation[atroposSeg == 6] <- 6
# Do left/right
if( verbose )
{
cat( " Doing left/right hemispheres.\n" )
}
leftLabels <- c( 4:7, 10:13, 17, 18, 25, 26, 28, 30, 91, 1002, 1003, 1005:1031, 1034, 1035 )
rightLabels <- c( 43:46, 49:54, 57, 58, 60, 62, 92, 2002, 2003, 2005:2031, 2034, 2035 )
hemisphereLabels <- list( leftLabels, rightLabels )
dktHemispheres <- antsImageClone( dktLabelImage )
for( i in seq.int( length( hemisphereLabels ) ) )
{
for( j in seq.int( length( hemisphereLabels[[i]] ) ) )
{
dktHemispheres[dktHemispheres == hemisphereLabels[[i]][j]] <- i
}
}
dktHemispheres[dktHemispheres > 2] <- 0
atroposBrainMask <- thresholdImage( atroposSeg, 0, 0, 0, 1 )
hemisphereParcellation <- iMath( atroposBrainMask, "PropagateLabelsThroughMask", atroposBrainMask * dktHemispheres )
for( i in seq.int( 6 ) )
{
lobarParcellation[lobarParcellation == i & hemisphereParcellation == 2] <- 6 + i
}
}
if( returnProbabilityImages && doLobarParcellation )
{
return( list(
segmentationImage = dktLabelImage,
lobarParcellation = lobarParcellation,
innerProbabilityImages = innerProbabilityImages,
outerProbabilityImages = outerProbabilityImages
)
)
} else if( returnProbabilityImages && ! doLobarParcellation ) {
return( list(
segmentationImage = dktLabelImage,
innerProbabilityImages = innerProbabilityImages,
outerProbabilityImages = outerProbabilityImages
)
)
} else if( ! returnProbabilityImages && doLobarParcellation ) {
return( list(
segmentationImage = dktLabelImage,
lobarParcellation = lobarParcellation
)
)
} else {
return( dktLabelImage )
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.