We build on the BRATS 2013 challenge to segment areas of the brain that have been damaged by stroke. We also refer to a more recent publication that implements a more complex version of what we do here.
We define a function that will help us simulate large, lateralized lesions on the fly.
library(ANTsR) simLesion<-function( img, s , w, thresh=0.01, mask=NA, myseed ) { set.seed(myseed) img<-iMath(img,"Normalize") if ( is.na(mask) ) mask<-getMask(img) i<-makeImage( dim(img) , rnorm( length(as.array(img)) ) ) i[ mask==0 ]<-0 ni<-smoothImage(i,s) ni[mask==0]<-0 i<-thresholdImage(ni,thresh,Inf) i<-iMath(i,"GetLargestComponent") ti<-antsImageClone(i) i[i>0]<-ti[i>0] i<-smoothImage(i,w) i[ mask != 1 ] <- 0 i[ 1:(dim(img)[1]/2), 1:(dim(img)[2]-1) ]<-0 limg<-( antsImageClone(img) * (-i) %>% iMath("Normalize") ) return( list(limg=limg, lesion=i ) ) }
Now let's apply this function to generate a test dataset.
ti<-antsImageRead( getANTsRData("r27") ) timask=getMask(ti) seg2<-kmeansSegmentation( ti, 3 )$segmentation ll2<-simLesion( ti, 10, 6, myseed=919 ) # different sized lesion seg2[ ll2$lesion > 0.5 & seg2 > 0.5 ]<-4
Now let's apply this function to generate a test dataset.
invisible( plot(ll2$limg) )
Create training data and map to the test subject. Note that
a "real" application of this type would use cost function masking.
But let's ignore that aspect of the problem here.
img<-antsImageRead( getANTsRData("r16") ) seg<-kmeansSegmentation( img, 3 )$segmentation ll<-simLesion( img, 12, 5, myseed=1 ) seg[ ll$lesion > 0.5 & seg > 0.5 ]<-4
invisible( plot(ll$limg) )
This gives us a subject with a "ground truth" segmentation.
Now we get a new subject and map to the space of the arbitrarily chosen reference space.
img<-antsImageRead( getANTsRData("r30") , 2 ) seg1<-kmeansSegmentation( img, 3 )$segmentation ll1<-simLesion( img, 9, 5, myseed=2 ) # different sized lesion seg1[ ll1$lesion > 0.5 & seg1 > 0.5 ]<-4
invisible( plot( ll1$limg ) )
Now use these to train a model.
rad<-c(1,1) # fast setting mr<-c(1,2,4,2,1) # multi-res schedule, U-style schedule masks=list( getMask(seg), getMask(seg1) ) rfm<-mrvnrfs( list(seg,seg1) , list(list(ll$limg), list(ll1$limg) ), masks, rad=rad, nsamples = 500, ntrees=1000, multiResSchedule=mr, voxchunk=500, do.trace = 100)
newrflist<-list() temp<-mrvnrfs( list(seg,seg1) , list(list(ll$limg), list(ll1$limg) ), masks, rad=rad, nsamples = 500, ntrees=1000, multiResSchedule=mr, voxchunk=500 ) for ( k in 1:length( mr ) ) if ( length( rfm$rflist[[k]]$classes ) == length( temp$rflist[[k]]$classes ) ) newrflist[[k]]<-combine( rfm$rflist[[k]], temp$rflist[[k]] ) rfm$rflist<-newrflist
We apply the learned model to segment the new data.
mmseg<-mrvnrfs.predict( rfm$rflist, list(list(ll2$limg)), timask, rad=rad, multiResSchedule=mr, voxchunk=500 )
invisible( plot(mmseg$seg[[1]], window.img=c(0,max(mmseg$seg[[1]]) ) ) )
Here is the ground truth.
invisible( plot( seg2, window.img=c(0,max(seg2) ) ) )
Take a quick look at the lesion probability.
invisible( plot( mmseg$probs[[1]][[ max(mmseg$seg[[1]]) ]] ) )
Now we compute the overlap.
dicenumer<-sum( mmseg$seg[[1]] == max(mmseg$seg[[1]]) & seg2 == max(seg2) ) dicedenom<-sum( mmseg$seg[[1]] == max(mmseg$seg[[1]]) ) + sum( seg2 == max(seg2) ) dice <- 2.0 * dicenumer / dicedenom if ( dice < 0.87 ) stop("suggests performance regression in mrvnrfs")
The Dice overlap is r dice
. We might consider model selection
as well where we do a quick estimate of lesion size based on the
volume of left hemisphere csf. Then build the model from
subjects that "match" with respect to the coarse amount of lesion.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.