diff --git a/NAMESPACE b/NAMESPACE index 74c7185..5da6158 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -124,7 +124,6 @@ export(decodeSsd3D) export(decodeUnet) export(deepAtropos) export(deepFlash) -export(deepFlashDeprecated) export(desikanKillianyTourvilleLabeling) export(drawRectangles) export(elBicho) diff --git a/R/deepFlash.R b/R/deepFlash.R index 4f0c1ab..99f3375 100644 --- a/R/deepFlash.R +++ b/R/deepFlash.R @@ -31,9 +31,11 @@ #' #' @param t1 raw or preprocessed 3-D T1-weighted brain image. #' @param t2 optional raw or preprocessed 3-D T2-weighted brain image. +#' @param whichParcellation string --- "yassa". See above label descriptions. #' @param doPreprocessing perform preprocessing. See description above. #' @param useRankIntensity If false, use histogram matching with cropped template -#' ROI. Otherwise, use a rank intensity transform on the cropped ROI. +#' ROI. Otherwise, use a rank intensity transform on the cropped ROI. Only for +#' 'yassa' parcellation. #' @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 @@ -51,846 +53,800 @@ #' results <- deepFlash( image ) #' } #' @export -deepFlash <- function( t1, t2 = NULL, doPreprocessing = TRUE, - antsxnetCacheDirectory = NULL, useRankIntensity = TRUE, verbose = FALSE ) +deepFlash <- function( t1, t2 = NULL, whichParcellation = "yassa", + doPreprocessing = TRUE, antsxnetCacheDirectory = NULL, + useRankIntensity = TRUE, verbose = FALSE ) { if( t1@dimension != 3 ) { stop( "Input image dimension must be 3." ) } - ################################ - # - # Options temporarily taken from the user - # - ################################ - - # @param useContralaterality If TRUE, use both hemispherical models to also - # predict the corresponding contralateral segmentation and use both sets of - # priors to produce the results. Mainly used for debugging. - # - # @param useHierarchicalParcellation If TRUE, use both hemispherical models to also - # predict the corresponding contralateral segmentation and use both sets of - # priors to produce the results. Mainly used for debugging. - - useHierarchicalParcellation <- TRUE - useContralaterality <- TRUE - - ################################ - # - # Preprocess image - # - ################################ - - t1Preprocessed <- t1 - t1Mask <- NULL - t1PreprocessedFlipped <- NULL - t1Template <- antsImageRead( getANTsXNetData( "deepFlashTemplateT1SkullStripped" ) ) - templateTransforms <- NULL - if( doPreprocessing ) + if( whichParcellation == "yassa" ) { - if( verbose ) - { - cat( "Preprocessing T1.\n" ) - } - - # Brain extraction - probabilityMask <- brainExtraction( t1Preprocessed, modality = "t1", - antsxnetCacheDirectory = antsxnetCacheDirectory, verbose = verbose ) - t1Mask <- thresholdImage( probabilityMask, 0.5, 1, 1, 0) - t1Preprocessed <- t1Preprocessed * t1Mask + ################################ + # + # Options temporarily taken from the user + # + ################################ - # Do bias correction - t1Preprocessed <- n4BiasFieldCorrection( t1Preprocessed, t1Mask, shrinkFactor = 4, verbose = verbose ) + # @param useContralaterality If TRUE, use both hemispherical models to also + # predict the corresponding contralateral segmentation and use both sets of + # priors to produce the results. Mainly used for debugging. + # + # @param useHierarchicalParcellation If TRUE, use both hemispherical models to also + # predict the corresponding contralateral segmentation and use both sets of + # priors to produce the results. Mainly used for debugging. - # Warp to template - registration <- antsRegistration( fixed = t1Template, moving = t1Preprocessed, - typeofTransform = "antsRegistrationSyNQuickRepro[a]", verbose = verbose ) - templateTransforms <- list( fwdtransforms = registration$fwdtransforms, - invtransforms = registration$invtransforms ) - t1Preprocessed <- registration$warpedmovout + useHierarchicalParcellation <- TRUE + useContralaterality <- TRUE - } - if( useContralaterality ) - { - t1PreprocessedDimension <- dim( t1Preprocessed ) - t1PreprocessedArray <- as.array( t1Preprocessed ) - t1PreprocessedArrayFlipped <- t1PreprocessedArray[t1PreprocessedDimension[1]:1,,] - t1PreprocessedFlipped <- as.antsImage( t1PreprocessedArrayFlipped, reference = t1Preprocessed ) - } + ################################ + # + # Preprocess image + # + ################################ - t2Preprocessed <- t2 - t2PreprocessedFlipped <- NULL - t2Template <- NULL - if( ! is.null( t2 ) ) - { - t2Template <- antsImageRead( getANTsXNetData( "deepFlashTemplateT2SkullStripped" ) ) - t2Template <- antsCopyImageInfo(t1Template, t2Template) + t1Preprocessed <- t1 + t1Mask <- NULL + t1PreprocessedFlipped <- NULL + t1Template <- antsImageRead( getANTsXNetData( "deepFlashTemplateT1SkullStripped" ) ) + templateTransforms <- NULL if( doPreprocessing ) { if( verbose ) { - cat( "Preprocessing T2.\n" ) + cat( "Preprocessing T1.\n" ) } # Brain extraction - t2Preprocessed <- t2Preprocessed * t1Mask + probabilityMask <- brainExtraction( t1Preprocessed, modality = "t1", + antsxnetCacheDirectory = antsxnetCacheDirectory, verbose = verbose ) + t1Mask <- thresholdImage( probabilityMask, 0.5, 1, 1, 0) + t1Preprocessed <- t1Preprocessed * t1Mask # Do bias correction - t2Preprocessed <- n4BiasFieldCorrection( t2Preprocessed, t1Mask, shrinkFactor = 4, verbose = verbose ) + t1Preprocessed <- n4BiasFieldCorrection( t1Preprocessed, t1Mask, shrinkFactor = 4, verbose = verbose ) # Warp to template - t2Preprocessed <- antsApplyTransforms( fixed = t1Template, moving = t2Preprocessed, - transformlist = templateTransforms$fwdtransforms, verbose = verbose ) + registration <- antsRegistration( fixed = t1Template, moving = t1Preprocessed, + typeofTransform = "antsRegistrationSyNQuickRepro[a]", verbose = verbose ) + templateTransforms <- list( fwdtransforms = registration$fwdtransforms, + invtransforms = registration$invtransforms ) + t1Preprocessed <- registration$warpedmovout + } if( useContralaterality ) { - t2PreprocessedDimension <- dim( t2Preprocessed ) - t2PreprocessedArray <- as.array( t2Preprocessed ) - t2PreprocessedArrayFlipped <- t2PreprocessedArray[t2PreprocessedDimension[1]:1,,] - t2PreprocessedFlipped <- as.antsImage(t2PreprocessedArrayFlipped, reference = t2Preprocessed ) + t1PreprocessedDimension <- dim( t1Preprocessed ) + t1PreprocessedArray <- as.array( t1Preprocessed ) + t1PreprocessedArrayFlipped <- t1PreprocessedArray[t1PreprocessedDimension[1]:1,,] + t1PreprocessedFlipped <- as.antsImage( t1PreprocessedArrayFlipped, reference = t1Preprocessed ) } - } - - probabilityImages <- list() - labels <- c( 0, 5:18 ) - imageSize <- c( 64, 64, 96 ) - - ################################ - # - # Process left/right in split network - # - ################################ - - ################################ - # - # Download spatial priors - # - ################################ - - spatialPriorsFileNamePath <- getANTsXNetData( "deepFlashPriors", - antsxnetCacheDirectory = antsxnetCacheDirectory ) - spatialPriors <- antsImageRead( spatialPriorsFileNamePath ) - priorsImageList <- splitNDImageToList( spatialPriors ) - for( i in seq.int( length( priorsImageList ) ) ) - { - priorsImageList[[i]] <- antsCopyImageInfo( t1Preprocessed, priorsImageList[[i]] ) - } - labelsLeft <- labels[seq.int( 2, length( labels ), by = 2)] - priorsImageLeftList <- priorsImageList[seq.int(2, length( priorsImageList ), by = 2 )] - probabilityImagesLeft <- list() - foregroundProbabilityImagesLeft <- list() - lowerBoundLeft <- c( 77, 75, 57 ) - upperBoundLeft <- c( 140, 138, 152 ) - tmpCropped <- cropIndices( t1Preprocessed, lowerBoundLeft, upperBoundLeft ) - originLeft <- antsGetOrigin( tmpCropped ) - - spacing <- antsGetSpacing( tmpCropped ) - direction <- antsGetDirection( tmpCropped ) - - t1TemplateRoiLeft <- cropIndices( t1Template, lowerBoundLeft, upperBoundLeft ) - t1TemplateRoiLeft <- iMath( t1TemplateRoiLeft, "Normalize" ) * 2.0 - 1.0 - t2TemplateRoiLeft <- NULL - if( ! is.null( t2Template ) ) - { - t2TemplateRoiLeft <- cropIndices( t2Template, lowerBoundLeft, upperBoundLeft ) - t2TemplateRoiLeft <- iMath( t2TemplateRoiLeft, "Normalize" ) * 2.0 - 1.0 - } - - labelsRight <- labels[seq.int( 3, length( labels ), by = 2)] - priorsImageRightList <- priorsImageList[seq.int(3, length( priorsImageList ), by = 2 )] - probabilityImagesRight <- list() - foregroundProbabilityImagesRight <- list() - lowerBoundRight <- c( 21, 75, 57 ) - upperBoundRight <- c( 84, 138, 152 ) - tmpCropped <- cropIndices( t1Preprocessed, lowerBoundRight, upperBoundRight ) - originRight <- antsGetOrigin( tmpCropped ) - - t1TemplateRoiRight <- cropIndices( t1Template, lowerBoundRight, upperBoundRight ) - t1TemplateRoiRight <- iMath( t1TemplateRoiRight, "Normalize" ) * 2.0 - 1.0 - t2TemplateRoiRight <- NULL - if( ! is.null( t2Template ) ) - { - t2TemplateRoiRight <- cropIndices( t2Template, lowerBoundRight, upperBoundRight ) - t2TemplateRoiRight <- iMath( t2TemplateRoiRight, "Normalize" ) * 2.0 - 1.0 - } + t2Preprocessed <- t2 + t2PreprocessedFlipped <- NULL + t2Template <- NULL + if( ! is.null( t2 ) ) + { + t2Template <- antsImageRead( getANTsXNetData( "deepFlashTemplateT2SkullStripped" ) ) + t2Template <- antsCopyImageInfo(t1Template, t2Template) + if( doPreprocessing ) + { + if( verbose ) + { + cat( "Preprocessing T2.\n" ) + } - ################################ - # - # Create model - # - ################################ + # Brain extraction + t2Preprocessed <- t2Preprocessed * t1Mask - channelSize <- 1 + length( labelsLeft ) - if( ! is.null( t2 ) ) - { - channelSize <- channelSize + 1 - } - numberOfClassificationLabels <- 1 + length( labelsLeft ) + # Do bias correction + t2Preprocessed <- n4BiasFieldCorrection( t2Preprocessed, t1Mask, shrinkFactor = 4, verbose = verbose ) - unetModel <- createUnetModel3D( c( imageSize, channelSize ), - numberOfOutputs = numberOfClassificationLabels, mode = 'classification', - numberOfFilters = c( 32, 64, 96, 128, 256 ), - convolutionKernelSize = c( 3, 3, 3 ), deconvolutionKernelSize = c( 2, 2, 2 ), - dropoutRate = 0.0, weightDecay = 0 ) + # Warp to template + t2Preprocessed <- antsApplyTransforms( fixed = t1Template, moving = t2Preprocessed, + transformlist = templateTransforms$fwdtransforms, verbose = verbose ) + } + if( useContralaterality ) + { + t2PreprocessedDimension <- dim( t2Preprocessed ) + t2PreprocessedArray <- as.array( t2Preprocessed ) + t2PreprocessedArrayFlipped <- t2PreprocessedArray[t2PreprocessedDimension[1]:1,,] + t2PreprocessedFlipped <- as.antsImage(t2PreprocessedArrayFlipped, reference = t2Preprocessed ) + } + } - penultimateLayer <- unetModel$layers[[length( unetModel$layers ) - 1]]$output + probabilityImages <- list() + labels <- c( 0, 5:18 ) + imageSize <- c( 64, 64, 96 ) - # medial temporal lobe - output1 <- penultimateLayer %>% keras::layer_conv_3d( filters = 1, kernel_size = 1L, - activation = 'sigmoid', kernel_regularizer = keras::regularizer_l2( 0.0 ) ) + ################################ + # + # Process left/right in split network + # + ################################ - if( useHierarchicalParcellation ) - { - # EC, perirhinal, and parahippo. - output2 <- penultimateLayer %>% keras::layer_conv_3d( filters = 1, kernel_size = 1L, - activation = 'sigmoid', kernel_regularizer = keras::regularizer_l2( 0.0 ) ) + ################################ + # + # Download spatial priors + # + ################################ - # Hippocampus - output3 <- penultimateLayer %>% keras::layer_conv_3d( filters = 1, kernel_size = 1L, - activation = 'sigmoid', kernel_regularizer = keras::regularizer_l2( 0.0 ) ) + spatialPriorsFileNamePath <- getANTsXNetData( "deepFlashPriors", + antsxnetCacheDirectory = antsxnetCacheDirectory ) + spatialPriors <- antsImageRead( spatialPriorsFileNamePath ) + priorsImageList <- splitNDImageToList( spatialPriors ) + for( i in seq.int( length( priorsImageList ) ) ) + { + priorsImageList[[i]] <- antsCopyImageInfo( t1Preprocessed, priorsImageList[[i]] ) + } - unetModel <- keras::keras_model( inputs = unetModel$input, outputs = list( unetModel$output, output1, output2, output3 ) ) - } else { - unetModel <- keras::keras_model( inputs = unetModel$input, outputs = list( unetModel$output, output1 ) ) - } + labelsLeft <- labels[seq.int( 2, length( labels ), by = 2)] + priorsImageLeftList <- priorsImageList[seq.int(2, length( priorsImageList ), by = 2 )] + probabilityImagesLeft <- list() + foregroundProbabilityImagesLeft <- list() + lowerBoundLeft <- c( 77, 75, 57 ) + upperBoundLeft <- c( 140, 138, 152 ) + tmpCropped <- cropIndices( t1Preprocessed, lowerBoundLeft, upperBoundLeft ) + originLeft <- antsGetOrigin( tmpCropped ) + + spacing <- antsGetSpacing( tmpCropped ) + direction <- antsGetDirection( tmpCropped ) + + t1TemplateRoiLeft <- cropIndices( t1Template, lowerBoundLeft, upperBoundLeft ) + t1TemplateRoiLeft <- iMath( t1TemplateRoiLeft, "Normalize" ) * 2.0 - 1.0 + t2TemplateRoiLeft <- NULL + if( ! is.null( t2Template ) ) + { + t2TemplateRoiLeft <- cropIndices( t2Template, lowerBoundLeft, upperBoundLeft ) + t2TemplateRoiLeft <- iMath( t2TemplateRoiLeft, "Normalize" ) * 2.0 - 1.0 + } - ################################ - # - # Left: build model and load weights - # - ################################ + labelsRight <- labels[seq.int( 3, length( labels ), by = 2)] + priorsImageRightList <- priorsImageList[seq.int(3, length( priorsImageList ), by = 2 )] + probabilityImagesRight <- list() + foregroundProbabilityImagesRight <- list() + lowerBoundRight <- c( 21, 75, 57 ) + upperBoundRight <- c( 84, 138, 152 ) + tmpCropped <- cropIndices( t1Preprocessed, lowerBoundRight, upperBoundRight ) + originRight <- antsGetOrigin( tmpCropped ) + + t1TemplateRoiRight <- cropIndices( t1Template, lowerBoundRight, upperBoundRight ) + t1TemplateRoiRight <- iMath( t1TemplateRoiRight, "Normalize" ) * 2.0 - 1.0 + t2TemplateRoiRight <- NULL + if( ! is.null( t2Template ) ) + { + t2TemplateRoiRight <- cropIndices( t2Template, lowerBoundRight, upperBoundRight ) + t2TemplateRoiRight <- iMath( t2TemplateRoiRight, "Normalize" ) * 2.0 - 1.0 + } - networkName <- 'deepFlashLeftT1' - if( ! is.null( t2 ) ) - { - networkName <- 'deepFlashLeftBoth' - } + ################################ + # + # Create model + # + ################################ - if( useHierarchicalParcellation ) - { - networkName <- paste0( networkName, "Hierarchical" ) - } + channelSize <- 1 + length( labelsLeft ) + if( ! is.null( t2 ) ) + { + channelSize <- channelSize + 1 + } + numberOfClassificationLabels <- 1 + length( labelsLeft ) - if( useRankIntensity ) - { - networkName <- paste0( networkName, "_ri" ) - } + unetModel <- createUnetModel3D( c( imageSize, channelSize ), + numberOfOutputs = numberOfClassificationLabels, mode = 'classification', + numberOfFilters = c( 32, 64, 96, 128, 256 ), + convolutionKernelSize = c( 3, 3, 3 ), deconvolutionKernelSize = c( 2, 2, 2 ), + dropoutRate = 0.0, weightDecay = 0 ) - if( verbose ) - { - cat( "DeepFlash: retrieving model weights (left).\n" ) - } - weightsFileName <- getPretrainedNetwork( networkName, antsxnetCacheDirectory = antsxnetCacheDirectory ) - load_model_weights_hdf5( unetModel, filepath = weightsFileName ) + penultimateLayer <- unetModel$layers[[length( unetModel$layers ) - 1]]$output - ################################ - # - # Left: do prediction and normalize to native space - # - ################################ + # medial temporal lobe + output1 <- penultimateLayer %>% keras::layer_conv_3d( filters = 1, kernel_size = 1L, + activation = 'sigmoid', kernel_regularizer = keras::regularizer_l2( 0.0 ) ) - if( verbose ) - { - cat( "Prediction (left).\n" ) - } + if( useHierarchicalParcellation ) + { + # EC, perirhinal, and parahippo. + output2 <- penultimateLayer %>% keras::layer_conv_3d( filters = 1, kernel_size = 1L, + activation = 'sigmoid', kernel_regularizer = keras::regularizer_l2( 0.0 ) ) - batchX <- NULL - if( useContralaterality ) - { - batchX <- array( data = 0, dim = c( 2, imageSize, channelSize ) ) - } else { - batchX <- array( data = 0, dim = c( 1, imageSize, channelSize ) ) - } + # Hippocampus + output3 <- penultimateLayer %>% keras::layer_conv_3d( filters = 1, kernel_size = 1L, + activation = 'sigmoid', kernel_regularizer = keras::regularizer_l2( 0.0 ) ) - t1Cropped <- cropIndices( t1Preprocessed, lowerBoundLeft, upperBoundLeft ) - if( useRankIntensity ) - { - t1Cropped <- rankIntensity( t1Cropped ) - } else { - t1Cropped <- histogramMatchImage( t1Cropped, t1TemplateRoiLeft, 255, 64, FALSE ) - } - batchX[1,,,,1] <- as.array( t1Cropped ) - if( useContralaterality ) - { - t1Cropped <- cropIndices( t1PreprocessedFlipped, lowerBoundLeft, upperBoundLeft ) - if( useRankIntensity ) - { - t1Cropped <- rankIntensity( t1Cropped ) + unetModel <- keras::keras_model( inputs = unetModel$input, outputs = list( unetModel$output, output1, output2, output3 ) ) } else { - t1Cropped <- histogramMatchImage( t1Cropped, t1TemplateRoiLeft, 255, 64, FALSE ) + unetModel <- keras::keras_model( inputs = unetModel$input, outputs = list( unetModel$output, output1 ) ) } - batchX[2,,,,1] <- as.array( t1Cropped ) - } - if( ! is.null( t2 ) ) - { - t2Cropped <- cropIndices( t2Preprocessed, lowerBoundLeft, upperBoundLeft ) - if( useRankIntensity ) + + ################################ + # + # Left: build model and load weights + # + ################################ + + networkName <- 'deepFlashLeftT1' + if( ! is.null( t2 ) ) { - t2Cropped <- rankIntensity( t2Cropped ) - } else { - t2Cropped <- histogramMatchImage( t2Cropped, t2TemplateRoiLeft, 255, 64, FALSE ) + networkName <- 'deepFlashLeftBoth' } - batchX[1,,,,2] <- as.array( t2Cropped ) - if( useContralaterality ) + + if( useHierarchicalParcellation ) { - t2Cropped <- cropIndices( t2PreprocessedFlipped, lowerBoundLeft, upperBoundLeft ) - if( useRankIntensity ) - { - t2Cropped <- rankIntensity( t2Cropped ) - } else { - t2Cropped <- histogramMatchImage( t2Cropped, t2TemplateRoiLeft, 255, 64, FALSE ) - } - batchX[2,,,,2] <- as.array( t2Cropped ) + networkName <- paste0( networkName, "Hierarchical" ) } - } - for( i in seq.int( length( priorsImageLeftList ) ) ) - { - croppedPrior <- cropIndices( priorsImageLeftList[[i]], lowerBoundLeft, upperBoundLeft ) - for( j in seq.int( dim( batchX )[1] ) ) + if( useRankIntensity ) { - batchX[j,,,,i + ( channelSize - length( labelsLeft ) )] <- as.array( croppedPrior ) + networkName <- paste0( networkName, "_ri" ) } - } - - predictedData <- unetModel %>% predict( batchX, verbose = verbose ) - probabilityImagesList <- decodeUnet( predictedData[[1]], t1Cropped ) - for( i in seq.int( 1 + length( labelsLeft ) ) ) - { - for( j in seq.int( dim( predictedData[[1]] )[1] ) ) + if( verbose ) { - probabilityImage <- probabilityImagesList[[j]][[i]] - if( i > 1 ) - { - probabilityImage <- decropImage( probabilityImage, t1Preprocessed * 0 ) - } else { - probabilityImage <- decropImage( probabilityImage, t1Preprocessed * 0 + 1 ) - } - - if( j == 2 ) # flipped - { - probabilityArray <- as.array( probabilityImage ) - probabilityArrayDimension <- dim( probabilityImage ) - probabilityArrayFlipped <- probabilityArray[probabilityArrayDimension[1]:1,,] - probabilityImage <- as.antsImage( probabilityArrayFlipped, reference = probabilityImage ) - } - - if( doPreprocessing == TRUE ) - { - probabilityImage <- antsApplyTransforms( fixed = t1, moving = probabilityImage, - transformlist = templateTransforms$invtransforms, - whichtoinvert = c( TRUE ), interpolator = "linear", verbose = verbose ) - } - - if( j == 1 ) # not flipped - { - probabilityImagesLeft <- append( probabilityImagesLeft, probabilityImage ) - } else { # flipped - probabilityImagesRight <- append( probabilityImagesRight, probabilityImage ) - } + cat( "DeepFlash: retrieving model weights (left).\n" ) } - } + weightsFileName <- getPretrainedNetwork( networkName, antsxnetCacheDirectory = antsxnetCacheDirectory ) + load_model_weights_hdf5( unetModel, filepath = weightsFileName ) - ################################ - # - # Left: do prediction of mtl, hippocampal, and ec regions and normalize to native space - # - ################################ + ################################ + # + # Left: do prediction and normalize to native space + # + ################################ - for( i in seq.int( 2, length( predictedData ) ) ) - { - probabilityImagesList <- decodeUnet( predictedData[[i]], t1Cropped ) - for( j in seq.int( dim( predictedData[[i]] )[1] ) ) + if( verbose ) { - probabilityImage <- probabilityImagesList[[j]][[1]] - probabilityImage <- decropImage( probabilityImage, t1Preprocessed * 0 ) - - if( j == 2 ) # flipped - { - probabilityArray <- as.array( probabilityImage ) - probabilityArrayDimension <- dim( probabilityImage ) - probabilityArrayFlipped <- probabilityArray[probabilityArrayDimension[1]:1,,] - probabilityImage <- as.antsImage( probabilityArrayFlipped, reference = probabilityImage ) - } - - if( doPreprocessing ) - { - probabilityImage <- antsApplyTransforms( fixed = t1, moving = probabilityImage, - transformlist = templateTransforms$invtransforms, - whichtoinvert = c( TRUE ), interpolator = "linear", verbose = verbose ) - } - - if( j == 1 ) # not flipped - { - foregroundProbabilityImagesLeft <- append( foregroundProbabilityImagesLeft, probabilityImage ) - } else { # flipped - foregroundProbabilityImagesRight <- append( foregroundProbabilityImagesRight, probabilityImage ) - } + cat( "Prediction (left).\n" ) } - } - - ################################ - # - # Right: build model and load weights - # - ################################ - networkName <- 'deepFlashRightT1' - if( ! is.null( t2 ) ) - { - networkName <- 'deepFlashRightBoth' - } - - if( useHierarchicalParcellation ) - { - networkName <- paste0( networkName, "Hierarchical" ) - } - - if( useRankIntensity ) - { - networkName <- paste0( networkName, "_ri" ) - } - - if( verbose ) - { - cat( "DeepFlash: retrieving model weights (Right).\n" ) - } - weightsFileName <- getPretrainedNetwork( networkName, antsxnetCacheDirectory = antsxnetCacheDirectory ) - load_model_weights_hdf5( unetModel, filepath = weightsFileName ) - - ################################ - # - # Right: do prediction and normalize to native space - # - ################################ - - if( verbose ) - { - cat( "Prediction (Right).\n" ) - } - - batchX <- NULL - if( useContralaterality ) - { - batchX <- array( data = 0, dim = c( 2, imageSize, channelSize ) ) - } else { - batchX <- array( data = 0, dim = c( 1, imageSize, channelSize ) ) - } - - t1Cropped <- cropIndices( t1Preprocessed, lowerBoundRight, upperBoundRight ) - if( useRankIntensity ) - { - t1Cropped <- rankIntensity( t1Cropped ) - } else { - t1Cropped <- histogramMatchImage( t1Cropped, t1TemplateRoiRight, 255, 64, FALSE ) - } - batchX[1,,,,1] <- as.array( t1Cropped ) - if( useContralaterality ) - { - t1Cropped <- cropIndices( t1PreprocessedFlipped, lowerBoundRight, upperBoundRight ) - if( useRankIntensity ) + batchX <- NULL + if( useContralaterality ) { - t1Cropped <- rankIntensity( t1Cropped ) + batchX <- array( data = 0, dim = c( 2, imageSize, channelSize ) ) } else { - t1Cropped <- histogramMatchImage( t1Cropped, t1TemplateRoiRight, 255, 64, FALSE ) + batchX <- array( data = 0, dim = c( 1, imageSize, channelSize ) ) } - batchX[2,,,,1] <- as.array( t1Cropped ) - } - if( ! is.null( t2 ) ) - { - t2Cropped <- cropIndices( t2Preprocessed, lowerBoundRight, upperBoundRight ) + + t1Cropped <- cropIndices( t1Preprocessed, lowerBoundLeft, upperBoundLeft ) if( useRankIntensity ) { - t2Cropped <- rankIntensity( t2Cropped ) + t1Cropped <- rankIntensity( t1Cropped ) } else { - t2Cropped <- histogramMatchImage( t2Cropped, t2TemplateRoiRight, 255, 64, FALSE ) + t1Cropped <- histogramMatchImage( t1Cropped, t1TemplateRoiLeft, 255, 64, FALSE ) } - batchX[1,,,,2] <- as.array( t2Cropped ) + batchX[1,,,,1] <- as.array( t1Cropped ) if( useContralaterality ) { - t2Cropped <- cropIndices( t2PreprocessedFlipped, lowerBoundRight, upperBoundRight ) + t1Cropped <- cropIndices( t1PreprocessedFlipped, lowerBoundLeft, upperBoundLeft ) if( useRankIntensity ) { - t2Cropped <- rankIntensity( t2Cropped ) + t1Cropped <- rankIntensity( t1Cropped ) } else { - t2Cropped <- histogramMatchImage( t2Cropped, t2TemplateRoiRight, 255, 64, FALSE ) + t1Cropped <- histogramMatchImage( t1Cropped, t1TemplateRoiLeft, 255, 64, FALSE ) } - batchX[2,,,,2] <- as.array( t2Cropped ) + batchX[2,,,,1] <- as.array( t1Cropped ) } - } - - for( i in seq.int( length( priorsImageRightList ) ) ) - { - croppedPrior <- cropIndices( priorsImageRightList[[i]], lowerBoundRight, upperBoundRight ) - for( j in seq.int( dim( batchX )[1] ) ) + if( ! is.null( t2 ) ) { - batchX[j,,,,i + ( channelSize - length( labelsRight ) )] <- as.array( croppedPrior ) - } - } - - predictedData <- unetModel %>% predict( batchX, verbose = verbose ) - probabilityImagesList <- decodeUnet( predictedData[[1]], t1Cropped ) - - for( i in seq.int( 1 + length( labelsRight ) ) ) - { - for( j in seq.int( dim( predictedData[[1]] )[1] ) ) - { - probabilityImage <- probabilityImagesList[[j]][[i]] - if( i > 1 ) + t2Cropped <- cropIndices( t2Preprocessed, lowerBoundLeft, upperBoundLeft ) + if( useRankIntensity ) { - probabilityImage <- decropImage( probabilityImage, t1Preprocessed * 0 ) + t2Cropped <- rankIntensity( t2Cropped ) } else { - probabilityImage <- decropImage( probabilityImage, t1Preprocessed * 0 + 1 ) + t2Cropped <- histogramMatchImage( t2Cropped, t2TemplateRoiLeft, 255, 64, FALSE ) } - - if( j == 2 ) # flipped + batchX[1,,,,2] <- as.array( t2Cropped ) + if( useContralaterality ) { - probabilityArray <- as.array( probabilityImage ) - probabilityArrayDimension <- dim( probabilityImage ) - probabilityArrayFlipped <- probabilityArray[probabilityArrayDimension[1]:1,,] - probabilityImage <- as.antsImage( probabilityArrayFlipped, reference = probabilityImage ) + t2Cropped <- cropIndices( t2PreprocessedFlipped, lowerBoundLeft, upperBoundLeft ) + if( useRankIntensity ) + { + t2Cropped <- rankIntensity( t2Cropped ) + } else { + t2Cropped <- histogramMatchImage( t2Cropped, t2TemplateRoiLeft, 255, 64, FALSE ) + } + batchX[2,,,,2] <- as.array( t2Cropped ) } + } - if( doPreprocessing ) + for( i in seq.int( length( priorsImageLeftList ) ) ) + { + croppedPrior <- cropIndices( priorsImageLeftList[[i]], lowerBoundLeft, upperBoundLeft ) + for( j in seq.int( dim( batchX )[1] ) ) { - probabilityImage <- antsApplyTransforms( fixed = t1, moving = probabilityImage, - transformlist = templateTransforms$invtransforms, - whichtoinvert = c( TRUE ), interpolator = "linear", verbose = verbose ) + batchX[j,,,,i + ( channelSize - length( labelsLeft ) )] <- as.array( croppedPrior ) } + } + + predictedData <- unetModel %>% predict( batchX, verbose = verbose ) + probabilityImagesList <- decodeUnet( predictedData[[1]], t1Cropped ) - if( j == 1 ) # not flipped + for( i in seq.int( 1 + length( labelsLeft ) ) ) + { + for( j in seq.int( dim( predictedData[[1]] )[1] ) ) { - if( useContralaterality ) + probabilityImage <- probabilityImagesList[[j]][[i]] + if( i > 1 ) { - probabilityImagesRight[[i]] <- ( probabilityImagesRight[[i]] + probabilityImage ) / 2 + probabilityImage <- decropImage( probabilityImage, t1Preprocessed * 0 ) } else { + probabilityImage <- decropImage( probabilityImage, t1Preprocessed * 0 + 1 ) + } + + if( j == 2 ) # flipped + { + probabilityArray <- as.array( probabilityImage ) + probabilityArrayDimension <- dim( probabilityImage ) + probabilityArrayFlipped <- probabilityArray[probabilityArrayDimension[1]:1,,] + probabilityImage <- as.antsImage( probabilityArrayFlipped, reference = probabilityImage ) + } + + if( doPreprocessing ) + { + probabilityImage <- antsApplyTransforms( fixed = t1, moving = probabilityImage, + transformlist = templateTransforms$invtransforms, + whichtoinvert = c( TRUE ), interpolator = "linear", verbose = verbose ) + } + + if( j == 1 ) # not flipped + { + probabilityImagesLeft <- append( probabilityImagesLeft, probabilityImage ) + } else { # flipped probabilityImagesRight <- append( probabilityImagesRight, probabilityImage ) } - } else { # flipped - probabilityImagesLeft[[i]] <- ( probabilityImagesLeft[[i]] + probabilityImage ) / 2 } } - } - ################################ - # - # Right: do prediction of mtl, hippocampal, and ec regions and normalize to native space - # - ################################ + ################################ + # + # Left: do prediction of mtl, hippocampal, and ec regions and normalize to native space + # + ################################ - for( i in seq.int( 2, length( predictedData ) ) ) - { - probabilityImagesList <- decodeUnet( predictedData[[i]], t1Cropped ) - for( j in seq.int( dim( predictedData[[i]] )[1] ) ) + for( i in seq.int( 2, length( predictedData ) ) ) { - probabilityImage <- probabilityImagesList[[j]][[1]] - probabilityImage <- decropImage( probabilityImage, t1Preprocessed * 0 ) - - if( j == 2 ) # flipped + probabilityImagesList <- decodeUnet( predictedData[[i]], t1Cropped ) + for( j in seq.int( dim( predictedData[[i]] )[1] ) ) { - probabilityArray <- as.array( probabilityImage ) - probabilityArrayDimension <- dim( probabilityImage ) - probabilityArrayFlipped <- probabilityArray[probabilityArrayDimension[1]:1,,] - probabilityImage <- as.antsImage( probabilityArrayFlipped, reference = probabilityImage ) - } + probabilityImage <- probabilityImagesList[[j]][[1]] + probabilityImage <- decropImage( probabilityImage, t1Preprocessed * 0 ) - if( doPreprocessing ) - { - probabilityImage <- antsApplyTransforms( fixed = t1, moving = probabilityImage, - transformlist = templateTransforms$invtransforms, - whichtoinvert = c( TRUE ), interpolator = "linear", verbose = verbose ) - } + if( j == 2 ) # flipped + { + probabilityArray <- as.array( probabilityImage ) + probabilityArrayDimension <- dim( probabilityImage ) + probabilityArrayFlipped <- probabilityArray[probabilityArrayDimension[1]:1,,] + probabilityImage <- as.antsImage( probabilityArrayFlipped, reference = probabilityImage ) + } - if( j == 1 ) # not flipped - { - if( useContralaterality ) + if( doPreprocessing ) { - foregroundProbabilityImagesRight[[i-1]] <- ( foregroundProbabilityImagesRight[[i-1]] + probabilityImage ) / 2 - } else { + probabilityImage <- antsApplyTransforms( fixed = t1, moving = probabilityImage, + transformlist = templateTransforms$invtransforms, + whichtoinvert = c( TRUE ), interpolator = "linear", verbose = verbose ) + } + + if( j == 1 ) # not flipped + { + foregroundProbabilityImagesLeft <- append( foregroundProbabilityImagesLeft, probabilityImage ) + } else { # flipped foregroundProbabilityImagesRight <- append( foregroundProbabilityImagesRight, probabilityImage ) } - } else { # flipped - foregroundProbabilityImagesLeft[[i-1]] <- ( foregroundProbabilityImagesLeft[[i-1]] + probabilityImage ) / 2 } } - } - ################################ - # - # Combine priors - # - ################################ + ################################ + # + # Right: build model and load weights + # + ################################ - probabilityBackgroundImage <- antsImageClone( t1 ) * 0 - for( i in seq.int( from = 2, to = length( probabilityImagesLeft ) ) ) - { - probabilityBackgroundImage <- probabilityBackgroundImage + probabilityImagesLeft[[i]] - } - for( i in seq.int( from = 2, to = length( probabilityImagesRight ) ) ) - { - probabilityBackgroundImage <- probabilityBackgroundImage + probabilityImagesRight[[i]] - } + networkName <- 'deepFlashRightT1' + if( ! is.null( t2 ) ) + { + networkName <- 'deepFlashRightBoth' + } - count <- 1 - probabilityImages[[count]] <- probabilityBackgroundImage * -1 + 1 - count <- count + 1 - for( i in seq.int( from = 2, to = length( probabilityImagesLeft ) ) ) - { - probabilityImages[[count]] <- probabilityImagesLeft[[i]] - count <- count + 1 - probabilityImages[[count]] <- probabilityImagesRight[[i]] - count <- count + 1 - } + if( useHierarchicalParcellation ) + { + networkName <- paste0( networkName, "Hierarchical" ) + } - ################################ - # - # Convert probability images to segmentation - # - ################################ - - imageMatrix <- imageListToMatrix( probabilityImages[2:length( probabilityImages )], t1 * 0 + 1 ) - backgroundForegroundMatrix <- rbind( imageListToMatrix( list( probabilityImages[[1]] ), t1 * 0 + 1 ), - colSums( imageMatrix ) ) - foregroundMatrix <- matrix( apply( backgroundForegroundMatrix, 2, which.max ), nrow = 1 ) - 1 - segmentationMatrix <- ( matrix( apply( imageMatrix, 2, which.max ), nrow = 1 ) + 1 ) * foregroundMatrix - segmentationImage <- matrixToImages( segmentationMatrix, t1 * 0 + 1 )[[1]] - - relabeledImage <- antsImageClone( segmentationImage ) - for( i in seq.int( length( labels ) ) ) - { - relabeledImage[( segmentationImage == i )] <- labels[i] - } + if( useRankIntensity ) + { + networkName <- paste0( networkName, "_ri" ) + } - foregroundProbabilityImages <- list() - for( i in seq.int( length( foregroundProbabilityImagesLeft ) ) ) - { - foregroundProbabilityImages[[i]] <- foregroundProbabilityImagesLeft[[i]] + foregroundProbabilityImagesRight[[i]] - } + if( verbose ) + { + cat( "DeepFlash: retrieving model weights (Right).\n" ) + } + weightsFileName <- getPretrainedNetwork( networkName, antsxnetCacheDirectory = antsxnetCacheDirectory ) + load_model_weights_hdf5( unetModel, filepath = weightsFileName ) - results <- NULL - if( useHierarchicalParcellation ) - { - results <- list( segmentationImage = relabeledImage, - probabilityImages = probabilityImages, - medialTemporalLobelProbabilityImage = foregroundProbabilityImages[[1]], - otherRegionProbabilityImage = foregroundProbabilityImages[[2]], - hippocampalProbabilityImage = foregroundProbabilityImages[[3]] - ) - } else { - results <- list( segmentationImage = relabeledImage, - probabilityImages = probabilityImages, - medialTemporalLobelProbabilityImage = foregroundProbabilityImages[[1]] - ) - } + ################################ + # + # Right: do prediction and normalize to native space + # + ################################ - return( results ) -} + if( verbose ) + { + cat( "Prediction (Right).\n" ) + } -#' Hippocampal/Enthorhinal segmentation using "Deep Flash" -#' -#' Perform hippocampal/entorhinal segmentation in T1 images using -#' labels from Mike Yassa's lab -#' -#' \url{https://faculty.sites.uci.edu/myassa/} -#' -#' The labeling is as follows: -#' \itemize{ -#' \item{Label 0 :}{background} -#' \item{Label 5 :}{left aLEC} -#' \item{Label 6 :}{right aLEC} -#' \item{Label 7 :}{left pMEC} -#' \item{Label 8 :}{right pMEC} -#' \item{Label 9 :}{left perirhinal} -#' \item{Label 10:}{right perirhinal} -#' \item{Label 11:}{left parahippocampal} -#' \item{Label 12:}{right parahippocampal} -#' \item{Label 13:}{left DG/CA3} -#' \item{Label 14:}{right DG/CA3} -#' \item{Label 15:}{left CA1} -#' \item{Label 16:}{right CA1} -#' \item{Label 17:}{left subiculum} -#' \item{Label 18:}{right subiculum} -#' } -#' -#' 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 doPerHemisphere If TRUE, do prediction based on separate networks per -#' hemisphere. Otherwise, use the single network trained for both hemispheres. -#' @param whichHemisphereModels "old" or "new". -#' @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 -#' inst/extdata/ subfolder of the ANTsRNet package. -#' @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 <- deepFlash( image ) -#' } -#' @export -deepFlashDeprecated <- function( t1, doPreprocessing = TRUE, doPerHemisphere = TRUE, - whichHemisphereModels = "new", antsxnetCacheDirectory = NULL, verbose = FALSE ) -{ - if( t1@dimension != 3 ) - { - stop( "Input image dimension must be 3." ) - } + batchX <- NULL + if( useContralaterality ) + { + batchX <- array( data = 0, dim = c( 2, imageSize, channelSize ) ) + } else { + batchX <- array( data = 0, dim = c( 1, imageSize, channelSize ) ) + } + + t1Cropped <- cropIndices( t1Preprocessed, lowerBoundRight, upperBoundRight ) + if( useRankIntensity ) + { + t1Cropped <- rankIntensity( t1Cropped ) + } else { + t1Cropped <- histogramMatchImage( t1Cropped, t1TemplateRoiRight, 255, 64, FALSE ) + } + batchX[1,,,,1] <- as.array( t1Cropped ) + if( useContralaterality ) + { + t1Cropped <- cropIndices( t1PreprocessedFlipped, lowerBoundRight, upperBoundRight ) + if( useRankIntensity ) + { + t1Cropped <- rankIntensity( t1Cropped ) + } else { + t1Cropped <- histogramMatchImage( t1Cropped, t1TemplateRoiRight, 255, 64, FALSE ) + } + batchX[2,,,,1] <- as.array( t1Cropped ) + } + if( ! is.null( t2 ) ) + { + t2Cropped <- cropIndices( t2Preprocessed, lowerBoundRight, upperBoundRight ) + if( useRankIntensity ) + { + t2Cropped <- rankIntensity( t2Cropped ) + } else { + t2Cropped <- histogramMatchImage( t2Cropped, t2TemplateRoiRight, 255, 64, FALSE ) + } + batchX[1,,,,2] <- as.array( t2Cropped ) + if( useContralaterality ) + { + t2Cropped <- cropIndices( t2PreprocessedFlipped, lowerBoundRight, upperBoundRight ) + if( useRankIntensity ) + { + t2Cropped <- rankIntensity( t2Cropped ) + } else { + t2Cropped <- histogramMatchImage( t2Cropped, t2TemplateRoiRight, 255, 64, FALSE ) + } + batchX[2,,,,2] <- as.array( t2Cropped ) + } + } - ################################ - # - # Preprocess image - # - ################################ + for( i in seq.int( length( priorsImageRightList ) ) ) + { + croppedPrior <- cropIndices( priorsImageRightList[[i]], lowerBoundRight, upperBoundRight ) + for( j in seq.int( dim( batchX )[1] ) ) + { + batchX[j,,,,i + ( channelSize - length( labelsRight ) )] <- as.array( croppedPrior ) + } + } - t1Preprocessed <- t1 - if( doPreprocessing == TRUE ) - { - 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 - } + predictedData <- unetModel %>% predict( batchX, verbose = verbose ) + probabilityImagesList <- decodeUnet( predictedData[[1]], t1Cropped ) - probabilityImages <- list() - labels <- c( 0, 5:18 ) + for( i in seq.int( 1 + length( labelsRight ) ) ) + { + for( j in seq.int( dim( predictedData[[1]] )[1] ) ) + { + probabilityImage <- probabilityImagesList[[j]][[i]] + if( i > 1 ) + { + probabilityImage <- decropImage( probabilityImage, t1Preprocessed * 0 ) + } else { + probabilityImage <- decropImage( probabilityImage, t1Preprocessed * 0 + 1 ) + } - ################################ - # - # Process left/right in same network - # - ################################ + if( j == 2 ) # flipped + { + probabilityArray <- as.array( probabilityImage ) + probabilityArrayDimension <- dim( probabilityImage ) + probabilityArrayFlipped <- probabilityArray[probabilityArrayDimension[1]:1,,] + probabilityImage <- as.antsImage( probabilityArrayFlipped, reference = probabilityImage ) + } - if( doPerHemisphere == FALSE ) - { + if( doPreprocessing ) + { + probabilityImage <- antsApplyTransforms( fixed = t1, moving = probabilityImage, + transformlist = templateTransforms$invtransforms, + whichtoinvert = c( TRUE ), interpolator = "linear", verbose = verbose ) + } + + if( j == 1 ) # not flipped + { + if( useContralaterality ) + { + probabilityImagesRight[[i]] <- ( probabilityImagesRight[[i]] + probabilityImage ) / 2 + } else { + probabilityImagesRight <- append( probabilityImagesRight, probabilityImage ) + } + } else { # flipped + probabilityImagesLeft[[i]] <- ( probabilityImagesLeft[[i]] + probabilityImage ) / 2 + } + } + } ################################ # - # Build model and load weights + # Right: do prediction of mtl, hippocampal, and ec regions and normalize to native space # ################################ - templateSize <- c( 160L, 192L, 160L ) + for( i in seq.int( 2, length( predictedData ) ) ) + { + probabilityImagesList <- decodeUnet( predictedData[[i]], t1Cropped ) + for( j in seq.int( dim( predictedData[[i]] )[1] ) ) + { + probabilityImage <- probabilityImagesList[[j]][[1]] + probabilityImage <- decropImage( probabilityImage, t1Preprocessed * 0 ) - 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( j == 2 ) # flipped + { + probabilityArray <- as.array( probabilityImage ) + probabilityArrayDimension <- dim( probabilityImage ) + probabilityArrayFlipped <- probabilityArray[probabilityArrayDimension[1]:1,,] + probabilityImage <- as.antsImage( probabilityArrayFlipped, reference = probabilityImage ) + } - if( verbose ) + if( doPreprocessing ) + { + probabilityImage <- antsApplyTransforms( fixed = t1, moving = probabilityImage, + transformlist = templateTransforms$invtransforms, + whichtoinvert = c( TRUE ), interpolator = "linear", verbose = verbose ) + } + + if( j == 1 ) # not flipped + { + if( useContralaterality ) + { + foregroundProbabilityImagesRight[[i-1]] <- ( foregroundProbabilityImagesRight[[i-1]] + probabilityImage ) / 2 + } else { + foregroundProbabilityImagesRight <- append( foregroundProbabilityImagesRight, probabilityImage ) + } + } else { # flipped + foregroundProbabilityImagesLeft[[i-1]] <- ( foregroundProbabilityImagesLeft[[i-1]] + probabilityImage ) / 2 + } + } + } + + ################################ + # + # Combine priors + # + ################################ + + probabilityBackgroundImage <- antsImageClone( t1 ) * 0 + for( i in seq.int( from = 2, to = length( probabilityImagesLeft ) ) ) { - cat( "DeepFlash: retrieving model weights.\n" ) + probabilityBackgroundImage <- probabilityBackgroundImage + probabilityImagesLeft[[i]] + } + for( i in seq.int( from = 2, to = length( probabilityImagesRight ) ) ) + { + probabilityBackgroundImage <- probabilityBackgroundImage + probabilityImagesRight[[i]] + } + + count <- 1 + probabilityImages[[count]] <- probabilityBackgroundImage * -1 + 1 + count <- count + 1 + for( i in seq.int( from = 2, to = length( probabilityImagesLeft ) ) ) + { + probabilityImages[[count]] <- probabilityImagesLeft[[i]] + count <- count + 1 + probabilityImages[[count]] <- probabilityImagesRight[[i]] + count <- count + 1 } - weightsFileName <- getPretrainedNetwork( "deepFlash", antsxnetCacheDirectory = antsxnetCacheDirectory ) - load_model_weights_hdf5( unetModel, filepath = weightsFileName ) ################################ # - # Do prediction and normalize to native space + # Convert probability images to segmentation # ################################ - if( verbose ) + imageMatrix <- imageListToMatrix( probabilityImages[2:length( probabilityImages )], t1 * 0 + 1 ) + backgroundForegroundMatrix <- rbind( imageListToMatrix( list( probabilityImages[[1]] ), t1 * 0 + 1 ), + colSums( imageMatrix ) ) + foregroundMatrix <- matrix( apply( backgroundForegroundMatrix, 2, which.max ), nrow = 1 ) - 1 + segmentationMatrix <- ( matrix( apply( imageMatrix, 2, which.max ), nrow = 1 ) + 1 ) * foregroundMatrix + segmentationImage <- matrixToImages( segmentationMatrix, t1 * 0 + 1 )[[1]] + + relabeledImage <- antsImageClone( segmentationImage ) + for( i in seq.int( length( labels ) ) ) { - cat( "Prediction.\n" ) + relabeledImage[( segmentationImage == i )] <- labels[i] } - croppedImage <- padOrCropImageToSize( t1Preprocessed, templateSize ) - imageArray <- as.array( croppedImage ) + foregroundProbabilityImages <- list() + for( i in seq.int( length( foregroundProbabilityImagesLeft ) ) ) + { + foregroundProbabilityImages[[i]] <- foregroundProbabilityImagesLeft[[i]] + foregroundProbabilityImagesRight[[i]] + } + + results <- NULL + if( useHierarchicalParcellation ) + { + results <- list( segmentationImage = relabeledImage, + probabilityImages = probabilityImages, + medialTemporalLobelProbabilityImage = foregroundProbabilityImages[[1]], + otherRegionProbabilityImage = foregroundProbabilityImages[[2]], + hippocampalProbabilityImage = foregroundProbabilityImages[[3]] + ) + } else { + results <- list( segmentationImage = relabeledImage, + probabilityImages = probabilityImages, + medialTemporalLobelProbabilityImage = foregroundProbabilityImages[[1]] + ) + } - batchX <- array( data = imageArray, dim = c( 1, templateSize, 1 ) ) - batchX <- ( batchX - mean( batchX ) ) / sd( batchX ) + return( results ) - predictedData <- unetModel %>% predict( batchX, verbose = verbose ) - probabilityImagesList <- decodeUnet( predictedData, croppedImage ) + } else if( whichParcellation == "wip" ) { + + useContralaterality <- TRUE - for( i in seq.int( length( probabilityImagesList[[1]] ) ) ) + ################################ + # + # Preprocess image + # + ################################ + + t1Preprocessed <- t1 + t1Mask <- NULL + t1PreprocessedFlipped <- NULL + t1Template <- antsImageRead( getANTsXNetData( "deepFlashTemplate2T1SkullStripped" ) ) + templateTransforms <- NULL + if( doPreprocessing ) { - if( i > 1 ) - { - decroppedImage <- decropImage( probabilityImagesList[[1]][[i]], t1Preprocessed * 0 ) - } else { - decroppedImage <- decropImage( probabilityImagesList[[1]][[i]], t1Preprocessed * 0 + 1 ) - } - if( doPreprocessing == TRUE ) + if( verbose ) { - probabilityImages[[i]] <- antsApplyTransforms( fixed = t1, moving = decroppedImage, - transformlist = t1Preprocessing$templateTransforms$invtransforms, - whichtoinvert = c( TRUE ), interpolator = "linear", verbose = verbose ) - } else { - probabilityImages[[i]] <- decroppedImage + cat( "Preprocessing T1.\n" ) } + + # Brain extraction + probabilityMask <- brainExtraction( t1Preprocessed, modality = "t1", + antsxnetCacheDirectory = antsxnetCacheDirectory, verbose = verbose ) + t1Mask <- thresholdImage( probabilityMask, 0.5, 1, 1, 0) + t1Preprocessed <- t1Preprocessed * t1Mask + + # Do bias correction + t1Preprocessed <- n4BiasFieldCorrection( t1Preprocessed, t1Mask, shrinkFactor = 4, verbose = verbose ) + + # Warp to template + registration <- antsRegistration( fixed = t1Template, moving = t1Preprocessed, + typeofTransform = "antsRegistrationSyNQuickRepro[a]", verbose = verbose ) + templateTransforms <- list( fwdtransforms = registration$fwdtransforms, + invtransforms = registration$invtransforms ) + t1Preprocessed <- registration$warpedmovout + + } + if( useContralaterality ) + { + t1PreprocessedDimension <- dim( t1Preprocessed ) + t1PreprocessedArray <- as.array( t1Preprocessed ) + t1PreprocessedArrayFlipped <- t1PreprocessedArray[t1PreprocessedDimension[1]:1,,] + t1PreprocessedFlipped <- as.antsImage( t1PreprocessedArrayFlipped, reference = t1Preprocessed ) } + probabilityImages <- list() + labelsLeft <- c( 104, 105, 106, 108, 109, 110, 114, 115, 126, 6001, 6003, 6008, 6009, 6010 ) + labelsRight <- c( 204, 205, 206, 208, 209, 210, 214, 215, 226, 7001, 7003, 7008, 7009, 7010 ) + + labels <- rep( 0, 1 + length( labelsLeft ) + length( labelsRight ) ) + labels[seq( 2, length( labels ), by = 2 )] <- labelsLeft + labels[seq( 3, length( labels ), by = 2 )] <- labelsRight + imageSize <- c( 64, 64, 128 ) + ################################ # # Process left/right in split networks # ################################ - } else { - ################################ # - # Left: download spatial priors + # Download spatial priors # ################################ - spatialPriorsLeftFileNamePath <- getANTsXNetData( "priorDeepFlashLeftLabels", + priorsLabelsFileNamePath <- getANTsXNetData( "deepFlashTemplate2Labels", antsxnetCacheDirectory = antsxnetCacheDirectory ) - spatialPriorsLeft <- antsImageRead( spatialPriorsLeftFileNamePath ) - priorsImageLeftList <- splitNDImageToList( spatialPriorsLeft ) + priorLabels <- antsImageRead( priorsLabelsFileNamePath ) + + priorsImageLeftList <- list() + for( i in seq.int( length( labelsLeft ) ) ) + { + priorImage <- thresholdImage( priorLabels, labelsLeft[i], labelsLeft[i], 1, 0 ) + priorsImage <- antsCopyImageInfo( t1Preprocessed, priorImage ) + priorsImageLeftList[[i]] <- smoothImage( priorImage, 1.0 ) + } + + priorsImageRightList <- list() + for( i in seq.int( length( labelsRight ) ) ) + { + priorImage <- thresholdImage( priorLabels, labelsRight[i], labelsRight[i], 1, 0 ) + priorsImage <- antsCopyImageInfo( t1Preprocessed, priorImage ) + priorsImageRightList[[i]] <- smoothImage( priorImage, 1.0 ) + } + + probabilityImagesLeft <- list() + foregroundProbabilityImagesLeft <- list() + lowerBoundLeft <- c( 115, 109, 83 ) + upperBoundLeft <- c( 178, 172, 210 ) + tmpCropped <- cropIndices( t1Preprocessed, lowerBoundLeft, upperBoundLeft ) + originLeft <- antsGetOrigin( tmpCropped ) + + spacing <- antsGetSpacing( tmpCropped ) + direction <- antsGetDirection( tmpCropped ) + + t1TemplateRoiLeft <- cropIndices( t1Template, lowerBoundLeft, upperBoundLeft ) + t1TemplateRoiLeft <- iMath( t1TemplateRoiLeft, "Normalize" ) * 2.0 - 1.0 + + probabilityImagesRight <- list() + foregroundProbabilityImagesRight <- list() + lowerBoundRight <- c( 51, 109, 83 ) + upperBoundRight <- c( 114, 172, 210 ) + tmpCropped <- cropIndices( t1Preprocessed, lowerBoundRight, upperBoundRight ) + originRight <- antsGetOrigin( tmpCropped ) + + t1TemplateRoiRight <- cropIndices( t1Template, lowerBoundRight, upperBoundRight ) + t1TemplateRoiRight <- iMath( t1TemplateRoiRight, "Normalize" ) * 2.0 - 1.0 ################################ # - # Left: build model and load weights + # Create model # ################################ - templateSize <- c( 64L, 96L, 96L ) - labelsLeft <- c( 0, 5, 7, 9, 11, 13, 15, 17 ) channelSize <- 1 + length( labelsLeft ) - - numberOfFilters <- 16 - networkName <- '' - if( whichHemisphereModels == "old" ) + if( ! is.null( t2 ) ) { - networkName <- "deepFlashLeft16" - } else if( whichHemisphereModels == "new" ) { - networkName <- "deepFlashLeft16new" - } else { - stop( "whichHemisphereModels must be \"old\" or \"new\"." ) + channelSize <- channelSize + 1 } + numberOfClassificationLabels <- 1 + length( labelsLeft ) - unetModel <- createUnetModel3D( c( templateSize, channelSize ), - numberOfOutputs = length( labelsLeft ), mode = 'classification', - numberOfLayers = 4, numberOfFiltersAtBaseLayer = numberOfFilters, dropoutRate = 0.0, + unetModel <- createUnetModel3D( c( imageSize, channelSize ), + numberOfOutputs = numberOfClassificationLabels, mode = 'classification', + numberOfFilters = c( 32, 64, 96, 128, 256 ), convolutionKernelSize = c( 3, 3, 3 ), deconvolutionKernelSize = c( 2, 2, 2 ), - weightDecay = 1e-5, additionalOptions = c( "attentionGating" ) ) + dropoutRate = 0.0, weightDecay = 0 ) + + penultimateLayer <- unetModel$layers[[length( unetModel$layers ) - 1]]$output + + # whole complex + output1 <- penultimateLayer %>% keras::layer_conv_3d( filters = 1, kernel_size = 1L, + activation = 'sigmoid', kernel_regularizer = keras::regularizer_l2( 0.0 ) ) + + # hippocampus + output2 <- penultimateLayer %>% keras::layer_conv_3d( filters = 1, kernel_size = 1L, + activation = 'sigmoid', kernel_regularizer = keras::regularizer_l2( 0.0 ) ) + + # amygdala + output3 <- penultimateLayer %>% keras::layer_conv_3d( filters = 1, kernel_size = 1L, + activation = 'sigmoid', kernel_regularizer = keras::regularizer_l2( 0.0 ) ) + + unetModel <- keras::keras_model( inputs = unetModel$input, + outputs = list( unetModel$output, output1, output2, output3 ) ) + + ################################ + # + # Left: build model and load weights + # + ################################ + + networkName <- 'deepFlash2LeftT1Hierarchical' if( verbose ) { @@ -910,51 +866,110 @@ deepFlashDeprecated <- function( t1, doPreprocessing = TRUE, doPerHemisphere = T cat( "Prediction (left).\n" ) } - croppedImage <- cropIndices( t1Preprocessed, c( 31, 52, 1 ), c( 94, 147, 96 ) ) - imageArray <- as.array( croppedImage ) - imageArray <- ( imageArray - mean( imageArray ) ) / sd( imageArray ) + batchX <- NULL + if( useContralaterality ) + { + batchX <- array( data = 0, dim = c( 2, imageSize, channelSize ) ) + } else { + batchX <- array( data = 0, dim = c( 1, imageSize, channelSize ) ) + } + + t1Cropped <- cropIndices( t1Preprocessed, lowerBoundLeft, upperBoundLeft ) + t1Cropped <- histogramMatchImage( t1Cropped, t1TemplateRoiLeft, 255, 64, FALSE ) - batchX <- array( data = 0, dim = c( 1, templateSize, channelSize ) ) - batchX[1,,,,1] <- imageArray + batchX[1,,,,1] <- as.array( t1Cropped ) + if( useContralaterality ) + { + t1Cropped <- cropIndices( t1PreprocessedFlipped, lowerBoundLeft, upperBoundLeft ) + t1Cropped <- histogramMatchImage( t1Cropped, t1TemplateRoiLeft, 255, 64, FALSE ) + batchX[2,,,,1] <- as.array( t1Cropped ) + } for( i in seq.int( length( priorsImageLeftList ) ) ) { - croppedPrior <- cropIndices( priorsImageLeftList[[i]], c( 31, 52, 1 ), c( 94, 147, 96 ) ) - batchX[1,,,,i+1] <- as.array( croppedPrior ) + croppedPrior <- cropIndices( priorsImageLeftList[[i]], lowerBoundLeft, upperBoundLeft ) + for( j in seq.int( dim( batchX )[1] ) ) + { + batchX[j,,,,i + ( channelSize - length( labelsLeft ) )] <- as.array( croppedPrior ) + } } predictedData <- unetModel %>% predict( batchX, verbose = verbose ) - probabilityImagesList <- decodeUnet( predictedData, croppedImage ) + probabilityImagesList <- decodeUnet( predictedData[[1]], t1Cropped ) - probabilityImagesLeft <- list() - for( i in seq.int( length( probabilityImagesList[[1]] ) ) ) + for( i in seq.int( 1 + length( labelsLeft ) ) ) { - if( i > 1 ) - { - decroppedImage <- decropImage( probabilityImagesList[[1]][[i]], t1Preprocessed * 0 ) - } else { - decroppedImage <- decropImage( probabilityImagesList[[1]][[i]], t1Preprocessed * 0 + 1 ) - } - if( doPreprocessing == TRUE ) + for( j in seq.int( dim( predictedData[[1]] )[1] ) ) { - probabilityImagesLeft[[i]] <- antsApplyTransforms( fixed = t1, moving = decroppedImage, - transformlist = t1Preprocessing$templateTransforms$invtransforms, - whichtoinvert = c( TRUE ), interpolator = "linear", verbose = verbose ) - } else { - probabilityImagesLeft[[i]] <- decroppedImage + probabilityImage <- probabilityImagesList[[j]][[i]] + if( i > 1 ) + { + probabilityImage <- decropImage( probabilityImage, t1Preprocessed * 0 ) + } else { + probabilityImage <- decropImage( probabilityImage, t1Preprocessed * 0 + 1 ) + } + + if( j == 2 ) # flipped + { + probabilityArray <- as.array( probabilityImage ) + probabilityArrayDimension <- dim( probabilityImage ) + probabilityArrayFlipped <- probabilityArray[probabilityArrayDimension[1]:1,,] + probabilityImage <- as.antsImage( probabilityArrayFlipped, reference = probabilityImage ) + } + + if( doPreprocessing ) + { + probabilityImage <- antsApplyTransforms( fixed = t1, moving = probabilityImage, + transformlist = templateTransforms$invtransforms, + whichtoinvert = c( TRUE ), interpolator = "linear", verbose = verbose ) + } + + if( j == 1 ) # not flipped + { + probabilityImagesLeft <- append( probabilityImagesLeft, probabilityImage ) + } else { # flipped + probabilityImagesRight <- append( probabilityImagesRight, probabilityImage ) + } } } ################################ # - # Right: download spatial priors + # Left: do prediction of whole, hippocampal, and amygdala regions and normalize to native space # ################################ - spatialPriorsRightFileNamePath <- getANTsXNetData( "priorDeepFlashRightLabels", - antsxnetCacheDirectory = antsxnetCacheDirectory ) - spatialPriorsRight <- antsImageRead( spatialPriorsRightFileNamePath ) - priorsImageRightList <- splitNDImageToList( spatialPriorsRight ) + for( i in seq.int( 2, length( predictedData ) ) ) + { + probabilityImagesList <- decodeUnet( predictedData[[i]], t1Cropped ) + for( j in seq.int( dim( predictedData[[i]] )[1] ) ) + { + probabilityImage <- probabilityImagesList[[j]][[1]] + probabilityImage <- decropImage( probabilityImage, t1Preprocessed * 0 ) + + if( j == 2 ) # flipped + { + probabilityArray <- as.array( probabilityImage ) + probabilityArrayDimension <- dim( probabilityImage ) + probabilityArrayFlipped <- probabilityArray[probabilityArrayDimension[1]:1,,] + probabilityImage <- as.antsImage( probabilityArrayFlipped, reference = probabilityImage ) + } + + if( doPreprocessing ) + { + probabilityImage <- antsApplyTransforms( fixed = t1, moving = probabilityImage, + transformlist = templateTransforms$invtransforms, + whichtoinvert = c( TRUE ), interpolator = "linear", verbose = verbose ) + } + + if( j == 1 ) # not flipped + { + foregroundProbabilityImagesLeft <- append( foregroundProbabilityImagesLeft, probabilityImage ) + } else { # flipped + foregroundProbabilityImagesRight <- append( foregroundProbabilityImagesRight, probabilityImage ) + } + } + } ################################ # @@ -962,26 +977,7 @@ deepFlashDeprecated <- function( t1, doPreprocessing = TRUE, doPerHemisphere = T # ################################ - templateSize <- c( 64L, 96L, 96L ) - labelsRight <- c( 0, 6, 8, 10, 12, 14, 16, 18 ) - channelSize <- 1 + length( labelsRight ) - - numberOfFilters <- 16 - networkName <- '' - if( whichHemisphereModels == "old" ) - { - networkName <- "deepFlashRight16" - } else if( whichHemisphereModels == "new" ) { - networkName <- "deepFlashRight16new" - } else { - stop( "whichHemisphereModels must be \"old\" or \"new\"." ) - } - - unetModel <- createUnetModel3D( c( templateSize, channelSize ), - numberOfOutputs = length( labelsRight ), mode = 'classification', - numberOfLayers = 4, numberOfFiltersAtBaseLayer = numberOfFilters, dropoutRate = 0.0, - convolutionKernelSize = c( 3, 3, 3 ), deconvolutionKernelSize = c( 2, 2, 2 ), - weightDecay = 1e-5, additionalOptions = c( "attentionGating" ) ) + networkName <- 'deepFlash2RightT1Hierarchical' if( verbose ) { @@ -1001,38 +997,117 @@ deepFlashDeprecated <- function( t1, doPreprocessing = TRUE, doPerHemisphere = T cat( "Prediction (Right).\n" ) } - croppedImage <- cropIndices( t1Preprocessed, c( 89, 52, 1 ), c( 152, 147, 96 ) ) - imageArray <- as.array( croppedImage ) - imageArray <- ( imageArray - mean( imageArray ) ) / sd( imageArray ) + batchX <- NULL + if( useContralaterality ) + { + batchX <- array( data = 0, dim = c( 2, imageSize, channelSize ) ) + } else { + batchX <- array( data = 0, dim = c( 1, imageSize, channelSize ) ) + } - batchX <- array( data = 0, dim = c( 1, templateSize, channelSize ) ) - batchX[1,,,,1] <- imageArray + t1Cropped <- cropIndices( t1Preprocessed, lowerBoundRight, upperBoundRight ) + t1Cropped <- histogramMatchImage( t1Cropped, t1TemplateRoiRight, 255, 64, FALSE ) + batchX[1,,,,1] <- as.array( t1Cropped ) + if( useContralaterality ) + { + t1Cropped <- cropIndices( t1PreprocessedFlipped, lowerBoundRight, upperBoundRight ) + t1Cropped <- histogramMatchImage( t1Cropped, t1TemplateRoiRight, 255, 64, FALSE ) + batchX[2,,,,1] <- as.array( t1Cropped ) + } for( i in seq.int( length( priorsImageRightList ) ) ) { - croppedPrior <- cropIndices( priorsImageRightList[[i]], c( 89, 52, 1 ), c( 152, 147, 96 ) ) - batchX[1,,,,i+1] <- as.array( croppedPrior ) + croppedPrior <- cropIndices( priorsImageRightList[[i]], lowerBoundRight, upperBoundRight ) + for( j in seq.int( dim( batchX )[1] ) ) + { + batchX[j,,,,i + ( channelSize - length( labelsRight ) )] <- as.array( croppedPrior ) + } } predictedData <- unetModel %>% predict( batchX, verbose = verbose ) - probabilityImagesList <- decodeUnet( predictedData, croppedImage ) + probabilityImagesList <- decodeUnet( predictedData[[1]], t1Cropped ) - probabilityImagesRight <- list() - for( i in seq.int( length( probabilityImagesList[[1]] ) ) ) + for( i in seq.int( 1 + length( labelsRight ) ) ) { - if( i > 1 ) + for( j in seq.int( dim( predictedData[[1]] )[1] ) ) { - decroppedImage <- decropImage( probabilityImagesList[[1]][[i]], t1Preprocessed * 0 ) - } else { - decroppedImage <- decropImage( probabilityImagesList[[1]][[i]], t1Preprocessed * 0 + 1 ) + probabilityImage <- probabilityImagesList[[j]][[i]] + if( i > 1 ) + { + probabilityImage <- decropImage( probabilityImage, t1Preprocessed * 0 ) + } else { + probabilityImage <- decropImage( probabilityImage, t1Preprocessed * 0 + 1 ) + } + + if( j == 2 ) # flipped + { + probabilityArray <- as.array( probabilityImage ) + probabilityArrayDimension <- dim( probabilityImage ) + probabilityArrayFlipped <- probabilityArray[probabilityArrayDimension[1]:1,,] + probabilityImage <- as.antsImage( probabilityArrayFlipped, reference = probabilityImage ) + } + + if( doPreprocessing ) + { + probabilityImage <- antsApplyTransforms( fixed = t1, moving = probabilityImage, + transformlist = templateTransforms$invtransforms, + whichtoinvert = c( TRUE ), interpolator = "linear", verbose = verbose ) + } + + if( j == 1 ) # not flipped + { + if( useContralaterality ) + { + probabilityImagesRight[[i]] <- ( probabilityImagesRight[[i]] + probabilityImage ) / 2 + } else { + probabilityImagesRight <- append( probabilityImagesRight, probabilityImage ) + } + } else { # flipped + probabilityImagesLeft[[i]] <- ( probabilityImagesLeft[[i]] + probabilityImage ) / 2 + } } - if( doPreprocessing == TRUE ) + } + + ################################ + # + # Right: do prediction of whole, hippocampal, and amygdala regions and normalize to native space + # + ################################ + + for( i in seq.int( 2, length( predictedData ) ) ) + { + probabilityImagesList <- decodeUnet( predictedData[[i]], t1Cropped ) + for( j in seq.int( dim( predictedData[[i]] )[1] ) ) { - probabilityImagesRight[[i]] <- antsApplyTransforms( fixed = t1, moving = decroppedImage, - transformlist = t1Preprocessing$templateTransforms$invtransforms, - whichtoinvert = c( TRUE ), interpolator = "linear", verbose = verbose ) - } else { - probabilityImagesRight[[i]] <- decroppedImage + probabilityImage <- probabilityImagesList[[j]][[1]] + probabilityImage <- decropImage( probabilityImage, t1Preprocessed * 0 ) + + if( j == 2 ) # flipped + { + probabilityArray <- as.array( probabilityImage ) + probabilityArrayDimension <- dim( probabilityImage ) + probabilityArrayFlipped <- probabilityArray[probabilityArrayDimension[1]:1,,] + probabilityImage <- as.antsImage( probabilityArrayFlipped, reference = probabilityImage ) + } + + if( doPreprocessing ) + { + probabilityImage <- antsApplyTransforms( fixed = t1, moving = probabilityImage, + transformlist = templateTransforms$invtransforms, + whichtoinvert = c( TRUE ), interpolator = "linear", verbose = verbose ) + } + + if( j == 1 ) # not flipped + { + if( useContralaterality ) + { + foregroundProbabilityImagesRight[[i-1]] <- ( foregroundProbabilityImagesRight[[i-1]] + probabilityImage ) / 2 + } else { + foregroundProbabilityImagesRight <- append( foregroundProbabilityImagesRight, probabilityImage ) + } + } else { # flipped + foregroundProbabilityImagesLeft[[i-1]] <- ( foregroundProbabilityImagesLeft[[i-1]] + probabilityImage ) / 2 + } } } @@ -1042,7 +1117,7 @@ deepFlashDeprecated <- function( t1, doPreprocessing = TRUE, doPerHemisphere = T # ################################ - probabilityBackgroundImage <- antsImageClone(t1) * 0 + probabilityBackgroundImage <- antsImageClone( t1 ) * 0 for( i in seq.int( from = 2, to = length( probabilityImagesLeft ) ) ) { probabilityBackgroundImage <- probabilityBackgroundImage + probabilityImagesLeft[[i]] @@ -1062,49 +1137,41 @@ deepFlashDeprecated <- function( t1, doPreprocessing = TRUE, doPerHemisphere = T probabilityImages[[count]] <- probabilityImagesRight[[i]] count <- count + 1 } - } - imageMatrix <- imageListToMatrix( probabilityImages[2:length( probabilityImages )], t1 * 0 + 1 ) - backgroundForegroundMatrix <- rbind( imageListToMatrix( list( probabilityImages[[1]] ), t1 * 0 + 1 ), - colSums( imageMatrix ) ) - foregroundMatrix <- matrix( apply( backgroundForegroundMatrix, 2, which.max ), nrow = 1 ) - 1 - segmentationMatrix <- ( matrix( apply( imageMatrix, 2, which.max ), nrow = 1 ) + 1 ) * foregroundMatrix - segmentationImage <- matrixToImages( segmentationMatrix, t1 * 0 + 1 )[[1]] + ################################ + # + # Convert probability images to segmentation + # + ################################ - relabeledImage <- antsImageClone( segmentationImage ) + imageMatrix <- imageListToMatrix( probabilityImages[2:length( probabilityImages )], t1 * 0 + 1 ) + backgroundForegroundMatrix <- rbind( imageListToMatrix( list( probabilityImages[[1]] ), t1 * 0 + 1 ), + colSums( imageMatrix ) ) + foregroundMatrix <- matrix( apply( backgroundForegroundMatrix, 2, which.max ), nrow = 1 ) - 1 + segmentationMatrix <- ( matrix( apply( imageMatrix, 2, which.max ), nrow = 1 ) + 1 ) * foregroundMatrix + segmentationImage <- matrixToImages( segmentationMatrix, t1 * 0 + 1 )[[1]] - for( i in seq.int( length( labels ) ) ) - { - relabeledImage[( segmentationImage == i )] <- labels[i] - } + relabeledImage <- antsImageClone( segmentationImage ) + for( i in seq.int( length( labels ) ) ) + { + relabeledImage[( segmentationImage == i )] <- labels[i] + } + + foregroundProbabilityImages <- list() + for( i in seq.int( length( foregroundProbabilityImagesLeft ) ) ) + { + foregroundProbabilityImages[[i]] <- foregroundProbabilityImagesLeft[[i]] + foregroundProbabilityImagesRight[[i]] + } - results <- list( segmentationImage = relabeledImage, probabilityImages = probabilityImages ) - - # debugging - - # if( debug == TRUE ) - # { - # inputImage <- unetModel$input - # featureLayer <- unetModel$layers[[length( unetModel$layers ) - 1]] - # featureFunction <- keras::backend()$`function`( list( inputImage ), list( featureLayer$output ) ) - # featureBatch <- featureFunction( list( batchX[1,,,,,drop = FALSE] ) ) - # featureImagesList <- decodeUnet( featureBatch[[1]], croppedImage ) - # featureImages <- list() - # for( i in seq.int( length( featureImagesList[[1]] ) ) ) - # { - # decroppedImage <- decropImage( featureImagesList[[1]][[i]], t1Preprocessed * 0 ) - # if( doPreprocessing == TRUE ) - # { - # featureImages[[i]] <- antsApplyTransforms( fixed = t1, moving = decroppedImage, - # transformlist = t1Preprocessing$templateTransforms$invtransforms, - # whichtoinvert = c( TRUE ), interpolator = "linear", verbose = verbose ) - # } else { - # featureImages[[i]] <- decroppedImage - # } - # } - # results[['featureImagesLastLayer']] <- featureImages - # } - - return( results ) + results <- list( segmentationImage = relabeledImage, + probabilityImages = probabilityImages, + wholeProbabilityImage = foregroundProbabilityImages[[1]], + hippocampalProbabilityImage = foregroundProbabilityImages[[2]], + amygdalaProbabilityImage = foregroundProbabilityImages[[3]] + ) + return( results ) + } else { + stop("Uncrecognized parcellation.") + } } diff --git a/R/getANTsXNetData.R b/R/getANTsXNetData.R index 5093a55..55c6527 100644 --- a/R/getANTsXNetData.R +++ b/R/getANTsXNetData.R @@ -29,6 +29,8 @@ getANTsXNetData <- function( "deepFlashTemplateT1SkullStripped", "deepFlashTemplateT2", "deepFlashTemplateT2SkullStripped", + "deepFlashTemplate2T1SkullStripped", + "deepFlashTemplate2Labels", "mprage_hippmapp3r", "protonLungTemplate", "ctLungTemplate", @@ -81,6 +83,8 @@ getANTsXNetData <- function( deepFlashTemplateT1SkullStripped = "https://figshare.com/ndownloader/files/31339867", deepFlashTemplateT2 = "https://figshare.com/ndownloader/files/31207798", deepFlashTemplateT2SkullStripped = "https://figshare.com/ndownloader/files/31339870", + deepFlashTemplate2T1SkullStripped = "https://figshare.com/ndownloader/files/46461451", + deepFlashTemplate2Labels = "https://figshare.com/ndownloader/files/46461415", mprage_hippmapp3r = "https://ndownloader.figshare.com/files/24984689", protonLungTemplate = "https://ndownloader.figshare.com/files/22707338", ctLungTemplate = "https://ndownloader.figshare.com/files/22707335", diff --git a/R/getPretrainedNetwork.R b/R/getPretrainedNetwork.R index 78b340e..9950c1c 100644 --- a/R/getPretrainedNetwork.R +++ b/R/getPretrainedNetwork.R @@ -73,6 +73,8 @@ getPretrainedNetwork <- function( "deepFlashRightT1Hierarchical_ri", "deepFlashLeftBothHierarchical_ri", "deepFlashRightBothHierarchical_ri", + "deepFlash2LeftT1Hierarchical", + "deepFlash2RightT1Hierarchical", "deepFlashLeft8", "deepFlashRight8", "deepFlashLeft16", @@ -226,6 +228,8 @@ getPretrainedNetwork <- function( deepFlashRightT1Hierarchical_ri = "https://figshare.com/ndownloader/files/33198800", deepFlashLeftBothHierarchical_ri = "https://figshare.com/ndownloader/files/33198803", deepFlashRightBothHierarchical_ri = "https://figshare.com/ndownloader/files/33198809", + deepFlash2LeftT1Hierarchical = "https://figshare.com/ndownloader/files/46461418", + deepFlash2RightT1Hierarchical = "https://figshare.com/ndownloader/files/46461421", denoising = "https://ndownloader.figshare.com/files/14235296", dktInner = "https://ndownloader.figshare.com/files/23266943", dktOuter = "https://ndownloader.figshare.com/files/23765132", diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml index 20bd8de..d0e6b3d 100644 --- a/docs/pkgdown.yml +++ b/docs/pkgdown.yml @@ -2,5 +2,5 @@ pandoc: 3.1.8 pkgdown: 2.0.7 pkgdown_sha: ~ articles: {} -last_built: 2024-08-12T22:11Z +last_built: 2024-08-23T21:32Z diff --git a/man/deepFlash.Rd b/man/deepFlash.Rd index 58dd208..dbf5843 100644 --- a/man/deepFlash.Rd +++ b/man/deepFlash.Rd @@ -7,6 +7,7 @@ deepFlash( t1, t2 = NULL, + whichParcellation = "yassa", doPreprocessing = TRUE, antsxnetCacheDirectory = NULL, useRankIntensity = TRUE, @@ -18,6 +19,8 @@ deepFlash( \item{t2}{optional raw or preprocessed 3-D T2-weighted brain image.} +\item{whichParcellation}{string --- "yassa". See above label descriptions.} + \item{doPreprocessing}{perform preprocessing. See description above.} \item{antsxnetCacheDirectory}{destination directory for storing the downloaded @@ -26,7 +29,8 @@ template and model weights. Since these can be resused, if inst/extdata/ subfolder of the ANTsRNet package.} \item{useRankIntensity}{If false, use histogram matching with cropped template -ROI. Otherwise, use a rank intensity transform on the cropped ROI.} +ROI. Otherwise, use a rank intensity transform on the cropped ROI. Only for +'yassa' parcellation.} \item{verbose}{print progress.} } diff --git a/man/deepFlashDeprecated.Rd b/man/deepFlashDeprecated.Rd deleted file mode 100644 index 214a643..0000000 --- a/man/deepFlashDeprecated.Rd +++ /dev/null @@ -1,85 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/deepFlash.R -\name{deepFlashDeprecated} -\alias{deepFlashDeprecated} -\title{Hippocampal/Enthorhinal segmentation using "Deep Flash"} -\usage{ -deepFlashDeprecated( - t1, - doPreprocessing = TRUE, - doPerHemisphere = TRUE, - whichHemisphereModels = "new", - antsxnetCacheDirectory = NULL, - verbose = FALSE -) -} -\arguments{ -\item{t1}{raw or preprocessed 3-D T1-weighted brain image.} - -\item{doPreprocessing}{perform preprocessing. See description above.} - -\item{doPerHemisphere}{If TRUE, do prediction based on separate networks per -hemisphere. Otherwise, use the single network trained for both hemispheres.} - -\item{whichHemisphereModels}{"old" or "new".} - -\item{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 -inst/extdata/ subfolder of the ANTsRNet package.} - -\item{verbose}{print progress.} -} -\value{ -list consisting of the segmentation image and probability images for -each label. -} -\description{ -Perform hippocampal/entorhinal segmentation in T1 images using -labels from Mike Yassa's lab -} -\details{ -\url{https://faculty.sites.uci.edu/myassa/} - -The labeling is as follows: -\itemize{ -\item{Label 0 :}{background} -\item{Label 5 :}{left aLEC} -\item{Label 6 :}{right aLEC} -\item{Label 7 :}{left pMEC} -\item{Label 8 :}{right pMEC} -\item{Label 9 :}{left perirhinal} -\item{Label 10:}{right perirhinal} -\item{Label 11:}{left parahippocampal} -\item{Label 12:}{right parahippocampal} -\item{Label 13:}{left DG/CA3} -\item{Label 14:}{right DG/CA3} -\item{Label 15:}{left CA1} -\item{Label 16:}{right CA1} -\item{Label 17:}{left subiculum} -\item{Label 18:}{right subiculum} -} - -Preprocessing on the training data consisted of: -\itemize{ -\item n4 bias correction, -\item denoising, -\item brain extraction, and -\item 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} -} -} -\examples{ -\dontrun{ -library( ANTsRNet ) -library( keras ) - -image <- antsImageRead( "t1.nii.gz" ) -results <- deepFlash( image ) -} -} -\author{ -Tustison NJ -}