Skull Stripping CT data
John Muschelli
2021-02-18
All code for this document is located at here.
Goal
In this tutorial, we will discuss skull-stripping (or brain-extracting) X-ray computed tomography (CT) scans. We will use data from TCIA (http://www.cancerimagingarchive.net/) as there is a great package called TCIApathfinder to interface with TCIA.
Using TCIApathfinder
In order to use TCIApathfinder, please see the vignette to obtain API keys. Here we will look at the collections:
library(TCIApathfinder)
library(dplyr)
collections = get_collection_names()
collections = collections$collection_names
head(collections)[1] "4D-Lung" "AAPM-RT-MAC" "ACRIN-DSC-MR-Brain" "ACRIN-FLT-Breast"
[5] "ACRIN-FMISO-Brain" "ACRIN-NSCLC-FDG-PET"mods = get_modality_names(body_part = "BREAST")
head(mods$modalities)[1] "CR" "CT" "MG" "MR" "OT" "PT"Getting Body Part Information
Here we can see all the parts of the body examined.
bp = get_body_part_names()
bp$body_parts [1] "ABD" "ABD PEL" "ABD PELV" "ABDOMEN"
[5] "ABDOMEN_PELVIS " "ABDOMENPELVIS" "AP PORTABLE CHE" "BD CT ABD WO_W "
[9] "BLADDER" "BRAIN" "BRAIN W/WO_AH32" "BREAST"
[13] "CAP" "CAROTID" "CERVIX" "CHEST"
[17] "CHEST (THORAX) " "CHEST COMPUTED " "CHEST NO GRID" "CHEST_ABDOMEN"
[21] "CHEST_TO_PELVIS" "CHEST/ABD" "CHESTABDOMEN" "CHESTABDPELVIS"
[25] "COLON" "CT 3PHASE REN" "CT CHEST W_ENHA" "CT CHEST WO CE"
[29] "CT THORAX W CNT" "CTA CHEST" "ESOPHAGUS" "EXTREMITY"
[33] "FUSION" "HEAD" "Head-and-Neck" "Head-Neck"
[37] "HEADANDNECK" "HEADNECK" "HEART" "J brzuszna"
[41] "J BRZUSZNA" "Kidney" "KIDNEY" "LEG"
[45] "LIVER" "LUMBO-SACRAL SP" "LUNG" "MEDIASTINUM"
[49] "NECK" "OUTSIDE FIL" "OVARY" "PANCREAS"
[53] "Pelvis" "PELVIS" "PET_ABDOMEN_PEL" "PET_CT SCAN CHE"
[57] "Phantom" "PHANTOM" "PORT CHEST" "PROSTATE"
[61] "RECTUM" "SEG" "SELLA" "SKULL"
[65] "SPI CHEST 5MM" "STOMACH" "TH CT CHEST WO " "Thorax"
[69] "THORAX" "THORAX CT _AH05" "THORAX CT _OT01" "THORAX_1HEAD_NE"
[73] "THORAXABD" "THYROID" "TSPINE" "UNDEFINED"
[77] "UTERUS" "WHOLEBODY" "WO INTER" Particularly, these areas are of interest. There seems to be a “bug” in TCIApathfinder::get_series_info which is acknowledged in the help file. Namely, the body_part_examined is not always a parameter to be set. We could get all the series info for all the collections from the code below, but it takes some times (> 15 minutes):
# could look for any of these
get_bp = c("BRAIN", "HEAD", "HEADNECK")
# takes a long time
res = pbapply::pblapply(collections, function(collection) {
x = get_series_info(
collection = collection,
modality = "CT")
x$series
})Getting Series
Here we will gather the series information for a study we know to have head CT data:
collection = "CPTAC-GBM"
series = get_series_info(
collection = collection,
modality = "CT")
series = series$series
head(series) patient_id collection study_instance_uid
1 NA CPTAC-GBM 1.3.6.1.4.1.14519.5.2.1.2857.3707.221249410799063035815783816913
2 NA CPTAC-GBM 1.3.6.1.4.1.14519.5.2.1.2857.3707.221249410799063035815783816913
3 NA CPTAC-GBM 1.3.6.1.4.1.14519.5.2.1.2857.3707.221249410799063035815783816913
4 NA CPTAC-GBM 1.3.6.1.4.1.14519.5.2.1.2857.3707.170705714007862724678123629040
5 NA CPTAC-GBM 1.3.6.1.4.1.14519.5.2.1.2857.3707.170705714007862724678123629040
6 NA CPTAC-GBM 1.3.6.1.4.1.14519.5.2.1.2857.3707.170705714007862724678123629040
series_instance_uid modality
1 1.3.6.1.4.1.14519.5.2.1.2857.3707.100565015879506080275493644685 CT
2 1.3.6.1.4.1.14519.5.2.1.2857.3707.176470763322052742670285487681 CT
3 1.3.6.1.4.1.14519.5.2.1.2857.3707.272098545527401893663335969793 CT
4 1.3.6.1.4.1.14519.5.2.1.2857.3707.254723691164851053423448594844 CT
5 1.3.6.1.4.1.14519.5.2.1.2857.3707.531177247834252562951224965872 CT
6 1.3.6.1.4.1.14519.5.2.1.2857.3707.225513954801691101397384975174 CT
protocol_name series_date series_description body_part_examined
1 1.6 CTA HEAD WITH WAND PROTOCOL 2001-01-15 SAG 10 X 2 MIP <NA>
2 1.6 CTA HEAD WITH WAND PROTOCOL 2001-01-15 AX 10 X 2 MIP <NA>
3 1.6 CTA HEAD WITH WAND PROTOCOL 2001-01-15 COR10 X 2 MIP <NA>
4 1.8 CTV HEAD Auto Transfer 75mL Iso 300 2001-01-23 CTV COR <NA>
5 1.8 CTV HEAD Auto Transfer 75mL Iso 300 2001-01-23 CTV SAG <NA>
6 1.8 CTV HEAD Auto Transfer 75mL Iso 300 2001-01-23 CTV AXIAL <NA>
series_number annotations_flag manufacturer manufacturer_model_name
1 603.000000 NA GE MEDICAL SYSTEMS LightSpeed VCT
2 601.000000 NA GE MEDICAL SYSTEMS LightSpeed VCT
3 602.000000 NA GE MEDICAL SYSTEMS LightSpeed VCT
4 602.000000 NA GE MEDICAL SYSTEMS LightSpeed VCT
5 603.000000 NA GE MEDICAL SYSTEMS LightSpeed VCT
6 601.000000 NA GE MEDICAL SYSTEMS LightSpeed VCT
software_versions image_count
1 <NA> 93
2 <NA> 124
3 <NA> 101
4 <NA> 107
5 <NA> 89
6 <NA> 53Here we grab the first series ID from this data which has a description of “HEAD STD” for standard head:
std_head = series %>%
filter(grepl("HEAD STD", series_description))
series_instance_uid = std_head$series_instance_uid[1]
download_unzip_series = function(series_instance_uid,
verbose = TRUE) {
tdir = tempfile()
dir.create(tdir, recursive = TRUE)
tfile = tempfile(fileext = ".zip")
tfile = basename(tfile)
if (verbose) {
message("Downloading Series")
}
res = save_image_series(
series_instance_uid = series_instance_uid,
out_dir = tdir,
out_file_name = tfile)
if (verbose) {
message("Unzipping Series")
}
stopifnot(file.exists(res$out_file))
tdir = tempfile()
dir.create(tdir, recursive = TRUE)
res = unzip(zipfile = res$out_file, exdir = tdir)
L = list(files = res,
dirs = unique(dirname(normalizePath(res))))
return(L)
}
# Download and unzip the image series
file_list = download_unzip_series(
series_instance_uid = series_instance_uid)Downloading SeriesUnzipping SeriesConverting DICOM to NIfTI
We will use dcm2niix to convert from DICOM to NIfTI. The function dcm2niix is wrapped in dcm2niir. We will use dcm2niir::dcm2nii to convert the file. We use check_dcm2nii to grab the relevant output files:
library(dcm2niir)
dcm_result = dcm2nii(file_list$dirs)#Copying Files# Converting to nii '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/dcm2niir/dcm2niix' -9 -v 1 -z y -f %p_%t_%s '/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpOK3Xb6/file5a1b6e441c7b'result = check_dcm2nii(dcm_result)Here we read the data into R into a nifti object:
library(neurobase)
img = readnii(result)
ortho2(img)range(img)[1] -3024 3071Here we will use neurobase::rescale_img to make sure the minimum is \(-1024\) and the maximum is \(3071\). The minimum can be lower for areas outside the field of view (FOV). Here we plot the image and the Winsorized version to see the brain tissue:
img = rescale_img(img, min.val = -1024, max.val = 3071)
ortho2(img)ortho2(img, window = c(0, 100))Skull Strip
We can skull strip the image using CT_Skull_Strip or CT_Skull_Stripper. The CT_Skull_Stripper has a simple switch to use CT_Skull_Strip or CT_Skull_Strip_robust.
library(ichseg)
ss = CT_Skull_Strip(img, verbose = FALSE)
ortho2(img, ss > 0,
window = c(0, 100),
col.y = scales::alpha("red", 0.5))The CT_Skull_Strip_robust function does 2 neck removals using remove_neck from extrantsr and then find the center of gravity (COG) twice to make sure the segmentation focuses on the head. In some instances, the whole neck is included in the scan, such as some of the head-neck studies in TCIA.
Showing a Robust Example with the neck
Getting Series
Here we will gather the series information for the Head-Neck Cetuximab collection:
collection = "Head-Neck Cetuximab"
series = get_series_info(
collection = collection,
modality = "CT")
series = series$series
whole_body = series %>%
filter(grepl("WB", series_description))file_list = download_unzip_series(
series_instance_uid = series$series_instance_uid[1])Downloading SeriesUnzipping Seriesdcm_result = dcm2nii(file_list$dirs, merge_files = TRUE)#Copying Files# Converting to nii '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/dcm2niir/dcm2niix' -9 -m y -v 1 -z y -f %p_%t_%s '/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpOK3Xb6/file5a1b27ae6d50'result = check_dcm2nii(dcm_result)Here we see the original data has a lot of the neck and some of the shoulders in the scan:
img = readnii(result)
img = rescale_img(img, min.val = -1024, max.val = 3071)
ortho2(img, window = c(0, 100))We will try CT_Skull_Strip without adding any robust options:
ss_wb = CT_Skull_Strip(img, verbose = FALSE)
ortho2(ss_wb, window = c(0, 100))We see that this does not work very well. We will use the robust version. Here we use CT_Skull_Stripper, which will call CT_Skull_Strip_robust. This will run extrantsr::remove_neck, runs CT_Skull_Strip, then estimates a new center of gravity (COG) and then run CT_Skull_Strip again, and then run some hole filling:
ss_wb_robust = CT_Skull_Stripper(img, verbose = FALSE, robust = TRUE)
ortho2(ss_wb_robust, window = c(0, 100))We see that this robust version works well for even data with the neck. We can try it on a whole body image as well.
The website data
We could also look at the website, but these do not always correspond to the API and get all the necessary results.
library(rvest)
library(dplyr)
x = read_html("https://www.cancerimagingarchive.net/collections/")
tab = html_table(x)[[1]]
head_tab = tab %>%
filter(grepl("Head|Brain", Location),
grepl("CT", `Image Types`),
Access == "Public")
brain_tab = tab %>%
filter(grepl("Brain", Location),
grepl("CT", `Image Types`),
Access == "Public")
brain_tab Collection Cancer Type Location Species
1 ACRIN-DSC-MR-Brain (ACRIN 6677/RTOG 0625) Glioblastoma Multiforme Brain Human
2 CPTAC-GBM Glioblastoma Multiforme Brain Human
3 ACRIN-FMISO-Brain (ACRIN 6684) Glioblastoma Brain Human
4 IvyGAP Glioblastoma Brain Human
5 TCGA-LGG Low Grade Glioma Brain Human
6 TCGA-GBM Glioblastoma Multiforme Brain Human
Subjects Image Types Supporting Data Access Status
1 123 MR, CT Clinical Public Complete
2 189 CT, CR, SC, MR, Pathology Clinical, Genomics, Proteomics Public Ongoing
3 45 CT, MR, PT Clinical Public Complete
4 39 MR, CT, Pathology Clinical, Genomics Public Complete
5 199 MR, CT, Pathology Clinical, Genomics, Image Analyses Public Complete
6 262 MR, CT, DX, Pathology Clinical, Genomics, Image Analyses Public Complete
Updated
1 2020-09-09
2 2020-03-31
3 2019-09-16
4 2016-12-30
5 2014-09-04
6 2014-05-08In brain_tab, we see we have a few collections. We are going to use the collection Head-Neck Cetuximab from above.
Getting Patient Information
We could sample patients from the collection here and get the patient information:
set.seed(20181203)
patients = get_patient_info(collection = collection)
info = patients$patients
head(info) patient_id patient_name patient_dob patient_sex patient_ethnic_group collection
1 0522c0001 0522c0001 NA F NA Head-Neck Cetuximab
2 0522c0002 0522c0002 NA M NA Head-Neck Cetuximab
3 0522c0003 0522c0003 NA M NA Head-Neck Cetuximab
4 0522c0009 0522c0009 NA M NA Head-Neck Cetuximab
5 0522c0013 0522c0013 NA M NA Head-Neck Cetuximab
6 0522c0070 0522c0070 NA M NA Head-Neck CetuximabThough we are not guaranteed the data will have Brain CT data. We will use the series variable to grab a relevant scan.