Simultaneous truth and performance level estimation (STAPLE) is a method for estimating a image segmentation. Commonly an image is segmented from a set of raters or readers. Here, we have an image that has been segmented from 3 different readers. Each is a binary segmentation, which includes areas of lesions from one person’s image.
Here we will use STAPLE to fuse the data:
library(stapler)
images = staple_example_data()
res = staple(x = images, set_orient = TRUE)
#> Reshaping images
#> All images are niftiImage ojects
#> Warning in `orientation<-`(`*tmp*`, value = ori): Image qform and sform codes
#> are both zero, so it cannot be reoriented
#> Warning in `orientation<-`(`*tmp*`, value = ori): Image qform and sform codes
#> are both zero, so it cannot be reoriented
#> Warning in `orientation<-`(`*tmp*`, value = ori): Image qform and sform codes
#> are both zero, so it cannot be reoriented
#> There are 2 levels present
#> iter: 1, diff: 0.881689818589628
#> Convergence!
#> Creating output probability images/arrays
#> Creating output label image/array
print(names(res))
#> [1] "sensitivity" "specificity" "probability"
#> [4] "label" "prior" "number_iterations"
#> [7] "convergence_threshold" "convergence_value" "converged"We use the set_orient = TRUE argument so that if the
images headers are different, then they will be set to the same
orientation when run. The staple function will read in the
images, reshape the data into a matrix, compute a prior for each element
(voxel, in this case), run STAPLE. The output probability image, labeled
image (if probability ≥ 0.5), and
prior image are given in the output res.
STAPLE works with multi-class data as well. The data do not need to be binary, but they need to be consistently labeled for each image/segmentation. Here we will
x = matrix(rbinom(5000, size = 5, prob = 0.5), ncol = 1000)
table(x)
#> x
#> 0 1 2 3 4 5
#> 162 787 1551 1545 798 157
res_mult = staple_multi_mat(x)
#> There are 6 levels present
#> Making multiple, matrices. Hot-one encode
#> iter: 25, diff: 1.46932572903102e-06
#> iter: 50, diff: 1.64688707116056e-10
#> iter: 75, diff: 2.00395255944841e-14
#> Convergence!
ncol(res_mult$probability)
#> [1] 6
colnames(res_mult$probability)
#> [1] "0" "1" "2" "3" "4" "5"