|
| 1 | +#' @title Whole-core frequency distribution of Hounsfield units |
| 2 | +#' |
| 3 | +#' @description Provides the raw data and plots a frequency distibution for Hounsfield Units in the entire core, also delineating material classes. |
| 4 | +#' |
| 5 | +#' @usage coreHist(directory = file.choose(), |
| 6 | +#' units = "percent", |
| 7 | +#' upperLim = 3045, lowerLim = -1025, |
| 8 | +#' airHU = -850.3233, airSD = 77.6953, |
| 9 | +#' SiHU = 271.7827, SiSD = 39.2814, |
| 10 | +#' glassHU = 1345.0696, glassSD = 45.4129, |
| 11 | +#' waterHU = 63.912, waterSD = 14.1728, |
| 12 | +#' returnData = TRUE, pngName = NULL) |
| 13 | +#' |
| 14 | +#' @param directory a character string that can be (1) a matrix of DICOM images that exists in the global environment, or (2) the address of an individual DICOM file in a folder of DICOM images. The default action is <code>file.choose()</code>; a browser menu appears so the user can select the the desired directory by identifying a single DICOM file in the folder of images. |
| 15 | +#' @param units units to be used for plotting purposes: either "percent" (the default) or "absolute" |
| 16 | +#' @param upperLim upper bound cutoff for pixels (Hounsfield Units); upper bound is inclusive |
| 17 | +#' @param lowerLim lower bound cutoff for pixels (Hounsfield Units); lower bound is exclusive |
| 18 | +#' @param airHU mean value for air-filled calibration rod (Hounsfield Units) |
| 19 | +#' @param airSD standard deviation for air-filled calibration rod |
| 20 | +#' @param SiHU mean value for colloidal silica calibration rod |
| 21 | +#' @param SiSD standard deviation for colloidal Si calibration rod |
| 22 | +#' @param glassHU mean value for glass calibration rod |
| 23 | +#' @param glassSD standard deviation for glass calibration rod |
| 24 | +#' @param waterHU mean value for water filled calibration rod |
| 25 | +#' @param waterSD standard deviation for water filled calibration rod |
| 26 | +#' @param returnData if \code{TRUE}, voxel counts for each Hounsfield unit from \code{lowerLim} to \code{upperLim} are returned, as are material class definitions. These are the data needed to re-create and modify the frequency plot. |
| 27 | +#' @param pngName if this is not \code{NULL}, the frequency plot is saved to disk. In that case, \code{pngName} should be a character string containing the name and address of the file. |
| 28 | +#' |
| 29 | +#' @return list if \code{returnData = TRUE}, a list is returned containing the frequencies for each Hounsfield unit value from \code{lowerLim} to \code{upperLim}, and (2) the boundaries for material classes. Lower boundaries for a component class are exclusive, while upper bounds are inclusive. These materials allow the frequency distribution to be plotted by the user. If \code{returnData = FALSE} the data are plotted in the graphics window, but nothing is preserved. |
| 30 | +#' |
| 31 | +#' @examples |
| 32 | +#' # data(core_426) |
| 33 | +#' coreHist("core_426", returnData = FALSE) |
| 34 | +#' |
| 35 | +#' @importFrom oro.dicom extractHeader |
| 36 | +#' @importFrom oro.dicom readDICOM |
| 37 | +#' @importFrom grDevices dev.off |
| 38 | +#' @importFrom grDevices png |
| 39 | +#' @importFrom graphics plot |
| 40 | +#' @importFrom graphics abline |
| 41 | +#' @importFrom graphics text |
| 42 | +#' @importFrom graphics par |
| 43 | +#' |
| 44 | +#' @export |
| 45 | + |
| 46 | +coreHist <- function(directory = file.choose(), |
| 47 | + units = "percent", |
| 48 | + upperLim = 3045, lowerLim = -1025, |
| 49 | + airHU = -850.3233, airSD = 77.6953, # all cal rod arguments are in Hounsfield Units |
| 50 | + SiHU = 271.7827, SiSD = 39.2814, |
| 51 | + glassHU = 1345.0696, glassSD = 45.4129, |
| 52 | + waterHU = 63.912, waterSD = 14.1728, |
| 53 | + returnData = TRUE, pngName = NULL) { |
| 54 | + if (!exists(directory)) { # is "directory" an existing object (user-loaded DICOM matrix) |
| 55 | + |
| 56 | + if (substr(directory, nchar(directory) - 3, nchar(directory)) %in% ".dcm") { |
| 57 | + directory <- dirname(directory) |
| 58 | + } else if (grep("/", directory) == 1) { # dangerously assumes that if there's a forward slash, it's a valid address |
| 59 | + # directory <- directory # do nothing |
| 60 | + } else stop("Incorrect directory name: directory specified in a character string must end with a '/'; if 'file.choose()' is used, the selected file must be a dicom image") |
| 61 | + # load DICOMs, takes a couple minutes |
| 62 | + fname <- oro.dicom::readDICOM(directory, verbose = TRUE) |
| 63 | + |
| 64 | + } else if (exists(directory) & (sum(names(get(directory)) %in% c("hdr", "img")) == 2)){ # could have better error checking here |
| 65 | + fname <- get(directory) |
| 66 | + } else stop("Invalid input: 'directory' object or file location is incorrectly specified.") |
| 67 | + |
| 68 | + # divisions between material classes |
| 69 | + splits <- data.frame(material = c("air", "RR", "water", "peat", "particles", "sand", "rock_shell"), |
| 70 | + lower = c(round(lowerLim), round(airHU + airSD), round(waterHU - waterSD), round(waterHU + waterSD), round(SiHU + SiSD), 750, round(glassHU + glassSD)), |
| 71 | + #lower = c(round(lowerLim), round(airHU+airSD) + 1, round(water.LB) + 1, round(water.UB) + 1, round(SiHU + SiSD) + 1, 750 + 1, round(glassHU + glassSD) + 1), |
| 72 | + upper = c(round(airHU + airSD), round(waterHU - waterSD), round(waterHU + waterSD), round(SiHU + SiSD), 750, round(glassHU + glassSD), round(upperLim))) |
| 73 | + |
| 74 | + |
| 75 | + pixelArea <- as.numeric(strsplit(fname$hdr[[1]]$value[fname$hdr[[1]]$name %in% "PixelSpacing"], " ")[[1]][1])^2 |
| 76 | + thick <- unique(oro.dicom::extractHeader(fname$hdr, "SliceThickness")) |
| 77 | + |
| 78 | + # convert raw units to Hounsfield units |
| 79 | + ct.slope <- unique(oro.dicom::extractHeader(fname$hdr, "RescaleSlope")) |
| 80 | + ct.int <- unique(oro.dicom::extractHeader(fname$hdr, "RescaleIntercept")) |
| 81 | + HU <- lapply(fname$img, function(x) x*ct.slope + ct.int) |
| 82 | + |
| 83 | + # tempDat <- as.data.frame(table(unlist(HU))) # Takes a long time. started at 7:27pm |
| 84 | + # tempDat <- lapply(HU, table) |
| 85 | + # tempDat2 <- plyr::join_all(lapply(tempDat, as.data.frame, stringsAsFactors = FALSE), by = "Var1") |
| 86 | + # tempDat2.temp <- do.call(rbind, HU) |
| 87 | + # tempDat2 <- data.frame(table(tempDat2.temp)) |
| 88 | + # names(tempDat2)[1] <- "Var1" |
| 89 | + # tempDat2$Var1 <- as.numeric(tempDat2$Var1) |
| 90 | + # tempDat2$finalFreq <- rowSums(tempDat2[, 2:ncol(tempDat2)], na.rm = TRUE) |
| 91 | + |
| 92 | + for ( i in 1:length(HU)) { |
| 93 | + if ( i == 1) { |
| 94 | + test <- tabulate((HU[[i]] - lowerLim), nbins = upperLim - lowerLim + 1) # all values need to be positive integers |
| 95 | + } else { |
| 96 | + test <- test + tabulate((HU[[i]] - lowerLim), nbins = upperLim - lowerLim + 1) |
| 97 | + } |
| 98 | + } |
| 99 | + tempDat2 <- data.frame(Var1 = (lowerLim + 1):(upperLim + 1), finalFreq = test) |
| 100 | + |
| 101 | + |
| 102 | + if (units == "percent") { |
| 103 | + tot <- sum(tempDat2$finalFreq, na.rm = TRUE) |
| 104 | + ylabel <- "Frequency (% of total voxels)" |
| 105 | + } else { |
| 106 | + tot <- 1 |
| 107 | + ylabel <- "Frequency (no. of voxels)" |
| 108 | + } |
| 109 | + if (!is.null(pngName)) { |
| 110 | + grDevices::png(file = pngName, width = 4, height = 3.5, units = "in") |
| 111 | + } |
| 112 | + graphics::par(mar = c(4, 4, 0.5, 0.5)) |
| 113 | + graphics::plot(finalFreq / tot ~ Var1, tempDat2[(tempDat2$Var1 < upperLim) & (tempDat2$Var1 > lowerLim), ], cex = 0.7, |
| 114 | + pch = 19, las = 1, ylab = ylabel, xlab = "HU", xlim = c(lowerLim, upperLim)) |
| 115 | + # add lines and label material classes |
| 116 | + graphics::abline(v = c(lowerLim + 1, splits$upper)) |
| 117 | + verticalTextPosition <- max(tempDat2[(tempDat2$Var1 < upperLim) & (tempDat2$Var1 > lowerLim), ], na.rm = TRUE) * c(0.90, 0.70) / tot |
| 118 | + graphics::text(x = (splits$upper + splits$lower)/2, y = verticalTextPosition, labels = as.character(splits$material)) # all classes |
| 119 | + # graphics::text(x = (splits$upper[-3] + splits$lower[-3])/2, y = verticalTextPosition, labels = as.character(splits$material)[-3]) # without water |
| 120 | + |
| 121 | + if (!is.null(pngName)) { |
| 122 | + grDevices::dev.off() |
| 123 | + } |
| 124 | + |
| 125 | + if (returnData == TRUE) { |
| 126 | + outDat <- tempDat2[(tempDat2$Var1 < upperLim) & (tempDat2$Var1 > lowerLim), c("Var1", "finalFreq")] |
| 127 | + return(list(histData = outDat, splits = splits)) |
| 128 | + } |
| 129 | + |
| 130 | +} |
0 commit comments