- Home
- Neuroconductor Tutorials
- Ss ct
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> 53
Here 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 Series
Unzipping Series
Converting 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 3071
Here 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))