Skull Stripping CT data

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))