Computed Tomography (CT) imaging has many applications outside the medical field. One such application in the field of environmental science is scanning sediment cores from coastal wetlands, a technique first demonstrated by Davey et al. (2011)[Davey, E., C. Wigand, R. Johnson, K. Sundberg, J. Morris, and C. Roman. 2011. Use of computed tomography imaging for quantifying coarse roots, rhizomes, peat, and particle densities in marsh soils. Ecological Applications 21(6): 2156-2171. (link to .pdf)]. The general goal is typically to quantify various soil/sediment components (particles, sand, water, roots) based on their densities; see sources citing Davey et al. (2011) for specific applications.
CT scanning produces large matrices of metadata-rich files with the .dcm extension. The extension is derived from data standard used in CT scanning, called Digital Imaging and Communications in Medicine (DICOM). The coreCT package is a small but powerful set of functions designed to streamline analysis of CT-scanned sediment cores, enabling rapid programmatic processing and synthesis of semi-processed DICOM images. 'Semi-processed' means that irrelevant parts of the images (PVC core tubes, etc.) have been masked out, and the only remaining data are for the actual sediment. coreCT builds on functionality in other packages, particulary the oro.dicom package (Whitcher et al. 2011)[Whitcher, B., V. J. Schmid and A. Thornton. 2011. Working with the DICOM and NIfTI Data Standards in R. Journal of Statistical Software 44(6): 1-28. (link)] for reading DICOM images and spatial analysis tools in the raster and igraph packages.
Key features of coreCT include:
This vignette introduces basic usage of the coreCT functions.
core_426 sample dataset included with coreCT includes three DICOM images from a CT-scanned sediment core. Each image represents a 0.625 mm depth interval, so this dataset is very small and only useful for demonstration purposes.
We'll use this provided dataset to demonstrate how the basic functions ,
rootSize work. The latter two functions work with a matrix of DICOM images already in the working environment, but coreCT also includes wrapper functions
rootSizeDir that are more flexible, operating on a directory of raw .dcm files and determining inputs like voxel volumes and slice thicknesses by automatically extracting metadata. These wrapper functions load the DICOM images, extract relevant metadata, and analyze the images for components and root characteristics.
Using the basic functions requires some initial work to extract metadata and convert the raw values to Hounsfield Units, the standard units for CT output.
This can be done more easily with coreCT::convDir, which automatically extracts DICOM metadata associated with pixel dimensions and converts raw values to Hounsfield Units:
materials <- convDir("core_426") plot(-depth ~ peat.cm3, data = materials, xlab = "Peat volume (cm3; per slice)", ylab = "Depth (cm)")
The conversion from Hounsfield Units to particle densities is done using calibration rods of known density analyzed with the sediment core. Reasonable default values are provided for four calibrants following Davey et al. (2011): air, water, colloidal silica, and glass. These values are used to define thresholds between soil components. These thresholds, and the distribution of values in the whole data series, can be visualized and returned by
HUfreq <- coreHist("core_426")
##  "histData" "splits"
## material lower upper ## 1 air -1025 -773 ## 2 RR -773 50 ## 3 water 50 78 ## 4 peat 78 311 ## 5 particles 311 750 ## 6 sand 750 1390 ## 7 rock_shell 1390 3045
conv uses these calibration data to classify voxels into material classes and then calculate four quantities for each of the seven material classes, for each image: (1) average Hounsfield units, (2) total 2D surface area (\(cm^2\)), (3) volume (\(cm^3\)), and (4) mass (g).
The 2D surface area (\(cm^2\)) is calculated as simply the number of pixels in each class, multiplied by the pixel area. Multiplying 2D surface area by the image thickness provides the volume of each component (reported in \(cm^3\)). The relationship between calibration rod Hounsfield units and particle densities is then applied to estimate component masses, by multiplying each voxel's bulk density by its volume (\(g.cm^3 * cm^3 = g\)).
Root characteristics are calculated by the
rootSize function, which also has a wrapper function
rootSizeDir that operates on a directory of raw DICOM images.
By default, root characteristics are calculated when
convDir is used, but if root data aren't of interest or a separate set of DICOM images is being used to quantify roots, time can be saved by including
rootData = FALSE as an argument to
Conversely, root characteristics can be analyzed without quantifying sediment composition. This uses the
rootSize function and its metadata-scraping wrapper function,
rootChars <- rootSizeDir("core_426", diameter.classes = c(1, 2.5, 10)) plot(-depth ~ structures, data = rootChars, xlab = "Root structures (per slice)", ylab = "Depth (cm)")