knitr::opts_chunk$set(echo = TRUE)

Load imaging data

DICOM files can be imported e.g. in MITK (http://www.mitk.org). MITK allows to segment whole brain and tumor volumes and then export of the resulting volumes as csv files.

These files usually look like this: x | y | z | VALUE whereas (x,y,z) are coordinates, and VALUE the respective value for T1ce, T2, .. sequences or CT data.

These files are imported in a first step. Several different variants are possible: - Import the raw csv files with usually 4 columns - Store (x,y,z) coordinates in one single value - +/- compression (e.g. using gzip)

The encoded coordinates are useful to easily perform set operations.

An encoded value has the following internal form, where {a,b,c} are the reserved bit positions for the {x,y,z} values.

x[x..x]ccccccccccbbbbbbbbbbaaaaaaaaaa

The number of repetitions is determined with the BITSIZE option. For the example data, 10 is used.

##TODO: check if this 'global var' can be avoided
options("BITSIZE"=10)

To ease analysis of imaging data, a new class is introducted: brainImage. This class inherits from the EBImage::Image class.

library(imageanalysisBrain)

For demonstration purposes, two brain datasets are included in this package, consisting of whole brain data an contoured tumor volumes.

The exemplary data is stored with an encoded coordinate (COORD).

First, both GTV and whole brain data is loaded.

## load GTV data, having only a COORD column (contains x,y,z) vals
gtvFile <- system.file("extdata", "TCGA-76-4931/PXL_GTV_COORD.csv.gz", 
    package = "imageanalysisBrain")
gtv <- read.table(gzfile(gtvFile), sep=",", header=TRUE)
gtvImage <- new("brainImage", gtv, gtvFile, 10) #10: Used Bitsize for encoding

##Load Brain data
brainFile <- system.file("extdata", "TCGA-76-4931/PXL_BRAIN_COORD.csv.gz", 
    package = "imageanalysisBrain")
brain <- read.table(gzfile(brainFile), sep=",", header=TRUE)
brainImage <- new("brainImage", brain, brainFile, 10) 

As different techniques are used for data acquisition, voxel sizes differ in their absolute dimensions. Here, an approximation is used which normalizes the maximum biparietal dimension of the whole brain data to 1. The scaling factor, which is later used to calculate distances, is stored in the xfactor slot.

## calculate x scaling factor 
brainImage <- imageanalysisBrain::calcXScalingFactor(brainImage)

Furthermore, no absolute center of origin within the whole brain volume is known. Therefore, a further heuristic is applied, which should localize a center of origin centrally between both lateral ventricles. This center of origin is automatically calculated and stored in the origin slot.

Image features

To explore characteristics of the available data, a wide range of image features can be extracted from the images.

These can be separated into different categories:

Technical differences in acquisitions between different datasets has to be compensated for.

The used heuristics are: - Tumor volumes are calculated as proportions of the whole brain volume - Distances are adjusted with the previously mentioned xfactor which normalized the maximal biparietal distance to 1.

## Store all imaging features
imgFeatures <- list()

# GTV in relation to whole Brain Volume
imgFeatures$gtvVolNorm <- 
    imageanalysisBrain::getNumberOfMeasurements(gtvImage)/
    imageanalysisBrain::getNumberOfMeasurements(brainImage)

# Only GTV dependent measurements
imgFeatures$gtvSD <- sd(imageanalysisBrain::getMeasurements(gtvImage))

#Distances between points of origins: GTv / Brain
imgFeatures$koordDistCanberraNorm <- dist(
    rbind(imageanalysisBrain::getOrigin(gtvImage)*
                imageanalysisBrain::getXFactor(brainImage), 
            imageanalysisBrain::getOrigin(brainImage)*
                imageanalysisBrain::getXFactor(brainImage)),
    method="canberra")[1]

Extract tissue classes from MRI data

To analyze tumor composition, brain / tumor tissue can be separated by intensity into different classes. However, absolute cutoffs, as can be used for CT images, are not applicable due to a often high hetereogeneity scanners / sequences.

Therefore, different "tissue classes" are extract by analyzing whole brain value to identify cutoffs for these classes. This allows to quantify the relative contibution of e.g. a hypointense class to a given tumor volume / GTV.

This separation is done by repeatedly drawing of small numbers of voxels from the whole brain volume. Their resultign distribution is then analyzed, and aggregated minimum values are selected as cutoffs.

Internally, a filtering step is applied, which selects only local maxima which are neighboured by two larger local maxima as valid cutoffs.

## takes approximately 40 sek for 100.000 iterations with 7 running threads
## on a i7-4820K CPU @ 3.70GHz GenuineIntel GNU/Linux
brainImage <- setThresholds(brainImage,
    imageanalysisBrain::identifyMinima(
        getMeasurements(brainImage), iterat=10000))

The complete dataset consisting of whole brain and tumor data can now be splitted into separated brainImage instances using the previously identified cutoffs. These can then be easily plotted (2D, 3D).

#TODO: Check for >= 1 subclasses ... 
subclassesGTV <- imageanalysisBrain::getTissueSubclasses(gtvImage, 
    imageanalysisBrain::getThresholds(brainImage))

subclassesBrain <- imageanalysisBrain::getTissueSubclasses(brainImage, 
    imageanalysisBrain::getThresholds(brainImage))

Plot brain / tumor slices

Single z-slices can be plotted as follows:

# Plot slice 10 (z) of the brain volume
EBImage::image(brainImage[,,10])
# and the GTV volume
EBImage::image(gtvImage[,,10])

If dimensions should be preserved to ensure direct comparability or overlap of images, a helper function (plotZSlice) can be used.

##Plot brain Slice
dim <- imageanalysisBrain::plotZSlice(brainImage, 10, dim=NULL) 

## ... and add the GTV slice 
imageanalysisBrain::plotZSlice(gtvImage, 10, 
    dim=dim, add=TRUE, col=gplots::bluered(1024)) 

Showing only subclass-slices works similar:

##Plot first brain subclass, z-slice 10
EBImage::image(subclassesBrain[[1]][,,10])


mknoll/imageanalysisBrain documentation built on May 23, 2019, 2:01 a.m.