GATE singles file processing

functions.R

filtering <- function(list, xFOV = 533, yFOV = 387, energy="All"){
  # unit of FOV is mm
  # energy windows including "Tc", "Ga", "I131"
  require(data.table)
  require(dplyr)
  require(magrittr)
  require(tidyr)
  require(Matrix)
  
  FOV <- subset(list, (V12 <= xFOV/2 & V12 >= -xFOV/2 & V14 <= yFOV/2 & V14 >= -yFOV/2))
  
  ###### energy window #####
  # unit of energy is MeV
  if (energy == "Tc"){
    print("You selected energy window of Tc-99m, 126.45~154.55 keV")
    filtered <- subset(FOV, V11>=0.12645 & V11 <=0.15455)
  }else if (energy == "Ga"){
    print("You selected energy window of Ga-67, 80.910 ~ 105.090, 165.6~202.4, 270~330 keV")
    filtered <- subset(FOV, (V11>=0.080910 & V11 <= 0.105090) |
                         (V11 >=0.1656 & V11<=0.2024)|
                         (V11 >=0.270 & V11 <= 0.330))
  }else if (energy == "I131"){
    print("You selected energy window of I-131, 327.6~400.4 keV")
    filtered <- subset(FOV, V11 >= 0.3276 & V11 <= 0.4004)
  }else{
    print("All energy are allowed")
    filtered <- FOV
  }
  
  return(filtered)
}

list2frame <- function(list, matrixSize=512, pixelSize = 1.12){
  # unit of pixel size is mm
  require(data.table)
  require(dplyr)
  require(magrittr)
  require(tidyr)
  require(Matrix)

  x <- list$V12 %/% pixelSize + matrixSize/2
  y <- list$V14 %/% pixelSize + matrixSize/2
  
  pixelArray <- data.frame(x,y) %>% count(x,y) %>%
    with(sparseMatrix(i=y,j=x, x=n, dims=c(matrixSize, matrixSize))) %>% ## notice the direction
    as.matrix
  
  image <- list(pixelSize = pixelSize, pixelArray=pixelArray)
  
  print(paste("pixel Size is", pixelSize, "; MatrixSize is", matrixSize))

  return(image)
}


energySpectrum <- function(list, EnergyLow = 0, EnergyHigh = 800, title){
  require(data.table)
  require(ggplot2)
  
  fig <- ggplot(list) + 
    geom_density(aes(x=V11*1000)) + 
    xlab("Energy (keV)") +
    xlim(EnergyLow, EnergyHigh) +
    labs(title=title)
  
  return(fig)
  
}

performance <- function(list){
  require(fitdistrplus)
  
  count <- dim(list)[1]
  
  fit_x <- fitdist(list$V12, "norm")
  fit_y <- fitdist(list$V14, "norm")
  
  FWHM_x <- fit_x$estimate[2]*2.355
  FWHM_y <- fit_y$estimate[2]*2.355
  
  result <- c(count, FWHM_x, FWHM_y)
  names(result) <-  c("count", "FWHM_x", "FWHM_y")
  
  return(result)
}