R/librMSIdata.R

In prafols/rMSI: MSI Tools

#########################################################################
#     rMSI - R package for MSI data handling and visualization
#     Copyright (C) 2014 Pere Rafols Soler
#
#     This program is free software: you can redistribute it and/or modify
#     it under the terms of the GNU General Public License as published by
#     the Free Software Foundation, either version 3 of the License, or
#     (at your option) any later version.
#
#     This program is distributed in the hope that it will be useful,
#     but WITHOUT ANY WARRANTY; without even the implied warranty of
#     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#     GNU General Public License for more details.
#
#     You should have received a copy of the GNU General Public License
#     along with this program.  If not, see <http://www.gnu.org/licenses/>.
############################################################################

#' Save a rMSI object to disk in a compressed .tar file.
#'
#' @param imgData a rMSI objecte with a image with a working ramdisk.
#' @param data_file Full path to hdd location to save image as a tar file.
#'
#' @export
SaveMsiData<-function(imgData, data_file)
{
  current_wd <- getwd() #Keep a copy of current working directory
  cat("Saving Image...\n")
  pt<-proc.time()

  setwd(dirname(data_file)) #Path of the resulting .tar file
  data_dir<-file.path(dirname(data_file), "ImgData")
  dir.create(data_dir)
  uuidObj<-imgData$uuid
  save(uuidObj, file = file.path(data_dir, "uuid.ImgR")) #Save UUID
  massObj<-imgData$mass
  save(massObj, file = file.path(data_dir, "mass.ImgR")) #Save mass axis
  sizeObj<-imgData$size
  save(sizeObj, file = file.path(data_dir, "size.ImgR")) #Save size Object
  posObj<-imgData$pos
  save(posObj, file = file.path(data_dir, "pos.ImgR")) #Save pos Object
  if(!is.null(imgData$posMotors))
  {
    posMotorsObj<-imgData$posMotors
    save(posMotorsObj, file = file.path(data_dir, "posMotors.ImgR")) #Save posMotors Object
  }
  meanSpcData<-imgData$mean
  save(meanSpcData, file = file.path(data_dir, "mean.SpcR")) #Save mean spectra
  resolutionObj<-imgData$pixel_size_um
  save(resolutionObj, file = file.path(data_dir, "pixel_size_um.ImgR")) #Save pixel size um Object
  if(!is.null(imgData$normalizations))
  {
    normalizationsObj <- imgData$normalizations
    save(normalizationsObj, file = file.path(data_dir, "normalizations.ImgR")) #Save normalizations Object
  }

  #Store also a vector of names of ff data in order to be able to restore it
  ffDataNames <- paste(names(imgData$data), "_ffzip", sep = "")
  save(ffDataNames, file = file.path(data_dir, "ffnames.ImgR")) #Save ff filenames Object

  ##New approach to avoid long and recursive path issues
  pb<-txtProgressBar(min = 0, max = length(ffDataNames), style = 3 )
  for(i in 1:length(ffDataNames))
  {
    setTxtProgressBar(pb, i)
    dm<-imgData$data[[i]][,]
    ffObj<-ff::ff(vmode = attr(attr(imgData$data[[i]],"physical"), "vmode"), dim = c(nrow(dm), ncol(dm)), filename =  file.path(data_dir, ffDataNames[i]))
    ffObj[,]<-dm
    ff::ffsave(ffObj , file =  file.path(data_dir, ffDataNames[i]))
    ff::delete(ffObj)
    rm(ffObj)
  }

  tar(tarfile =  file.path(getwd(), basename(data_file)), files = "ImgData")
  unlink(data_dir, recursive = T) #Remove intermediate data

  close(pb)
  pt<-proc.time() - pt
  cat(paste("Saving time:",round(pt["elapsed"], digits = 1),"seconds\n"))
  setwd(current_wd) #Restoire the original working directory
}


#' Load rMSI data from a compressed tar.
#'
#' @param data_file The tar o imzML file containing the MS image in rMSI format or imzML.
#' @param restore_path Where the ramdisk will be created.
#' @param fun_progress This is a callback function to update the progress of loading data. See details for more information.
#' @param ff_overwrite Tell ff to overwrite or not current ramdisk files.
#' @param fun_label This is a callback function to update the progress bar dialog text.
#' @param close_signal function to be called if the loading process is aborted.
#' @param imzMLChecksum if the binary file checksum must be verified, it can be disabled for convenice with really big files.
#' @param imzMLRename the image name, if NULL a default name based on the file name will be used.
#' @param imzMLSubCoords a Complex vector with the motors coordinates to be included in the ramdisk, if NULL all positions will be used.
#'
#' @return an rMSI object pointing to ramdisk stored data
#'
#' Loads a rMSI data object from .tar compressed file or imzML format. It will be uncompressed at specified restore_path.
#' fun_progress can be NULL or a function with the following prototipe: fun_progress( currentState ). If NULL is used
#' a default command line progress bar is used.
#' This function will be called periodically to monitor the loading status. This is usefull to implement progressbars.
#' If ramdisk is already created befor calling this method the parameter ff_overwrite will control the loadin behaviour. If it is set to false (default)
#' The ramdisk will be kept and the imaged loaded imediatelly. Otherwise if is set to true, the while dataset will be reloaded from tar file.
#'
#' @export
LoadMsiData<-function(data_file, restore_path = file.path(dirname(data_file), paste("ramdisk",basename(data_file), sep = "_")) , fun_progress = NULL, ff_overwrite = F, fun_label = NULL, close_signal = NULL, imzMLChecksum = F, imzMLRename = NULL, imzMLSubCoords = NULL)
{
  cat("Loading Image...\n")
  pt<-proc.time()

  #1- Check if the specified image ramdisk exists in the restore_path location
  datacube<-.FastLoad(restore_path)
  if(!is.null(datacube) && !ff_overwrite)
  {
    if(!is.null(fun_progress))
    {
      fun_progress(100)
    }

    pt<-proc.time() - pt
    cat(paste("Importing time:",round(pt["elapsed"], digits = 1),"seconds\n"))
    class(datacube) <- "rMSIObj"
    return(datacube)
  }

  fileExtension <- unlist(strsplit(basename(data_file), split = "\\."))
  fileExtension <- as.character(fileExtension[length(fileExtension)])
  if( fileExtension == "imzML")
  {
    return(import_imzML(data_file, ramdisk_path = restore_path, fun_progress = fun_progress, fun_text = fun_label, close_signal = close_signal, verifyChecksum = imzMLChecksum, subImg_rename = imzMLRename, subImg_Coords = imzMLSubCoords))
  }
  else if(fileExtension == "tar")
  {
    return(import_rMSItar(data_file,restore_path, fun_progress, fun_text = fun_label, close_signal = close_signal))
  }
  else
  {
    cat("The slected file is not valid.\n")
  }
  pt<-proc.time() - pt
  cat(paste("Importing time:",round(pt["elapsed"], digits = 1),"seconds\n"))
}

#' import_rMSItar.
#'
#' @param data_file The tar file containing the MS image in rMSI format.
#' @param restore_path Where the ramdisk will be created.
#' @param fun_progress This is a callback function to update the progress of loading data. See details for more information.
#' @param fun_text This is a callback function to update the label widget of loading data. See details for more information.
#' @param close_signal function to call if error.
#'
#'  Imports an rMSI data object in a tar data file
#'  It is recomanded to use rMSI::LoadMsiData directly instead of this function.
#'
#' @return  an rMSI data object.
#' @export
#'
import_rMSItar<-function(data_file, restore_path, fun_progress = NULL, fun_text = NULL, close_signal = NULL )
{
  setPbarValue<-function(progress)
  {
    setTxtProgressBar(pb, progress)
    return(T)
  }

  if(is.null(fun_progress))
  {
    pb<-txtProgressBar(min = 0, max = 100, style = 3 )
    fun_progress <- setPbarValue
  }
  else
  {
    pb<-NULL
  }

  #2 - Image is not preloaded so... the slow way
  if(!is.null(fun_text))
  {
    fun_text("Unpacking data...")
  }

  untar(tarfile = data_file, exdir = restore_path)
  img_path<-file.path(restore_path, "ImgData")


  if( file.exists(file.path(img_path, "uuid.ImgR") ))
  {
    load(file.path(img_path, "uuid.ImgR"))
  }
  else
  {
    #For compatibility with old datasets, just set no id
    uuidObj <- ""
  }
  
  load(file.path(img_path, "mass.ImgR"))
  load(file.path(img_path, "size.ImgR"))
  load(file.path(img_path, "pos.ImgR"))
  if( file.exists(file.path(img_path, "posMotors.ImgR")))
  {
    load(file.path(img_path, "posMotors.ImgR"))
  }
  else
  {
    posMotorsObj <- NULL
  }
  load(file.path(img_path, "mean.SpcR"))
  load(file.path(img_path, "ffnames.ImgR"))
  if(file.exists(file.path(img_path, "pixel_size_um.ImgR")))
  {
    load(file.path(img_path, "pixel_size_um.ImgR"))
  }
  else
  {
    print("Warning: Old image without resolution object. It is set to 9999 um by default!")
    resolutionObj <- 9999
  }
  if(file.exists(file.path(img_path, "normalizations.ImgR")))
  {
    load(file.path(img_path, "normalizations.ImgR"))
  }
  else
  {
    normalizationsObj <- NULL
  }

  if(!is.null(fun_text))
  {
    fun_text("Loading data in the ramdisk...")
  }
  spectra<-list()
  ppStep<-100/length(ffDataNames)
  pp<-0
  for(i in 1:length(ffDataNames))
  {
    ff::ffload(file.path(img_path, ffDataNames[i]), rootpath = restore_path, overwrite = T)
     #Get Hdd space back asap
    unlink(file.path(img_path, paste0(ffDataNames[i], ".RData")))
    unlink(file.path(img_path, paste0(ffDataNames[i], ".ffData")))
    spectra[[i]]<-ffObj
    names(spectra)[i]<-paste("ramdisk",i,".dat",sep = "")

    pp_ant<-pp
    pp<-pp+ppStep
    if(!is.null(fun_progress) && (round(pp) > round(pp_ant)) )
    {
      #Update progress bar
      if( !fun_progress(pp) )
      {
        return(NULL) #progress bar function must return true if the loading process is to be continued.
      }
    }
  }

  lapply(spectra, function(x){open(x)})
  datacube<-list(name = basename(data_file), uuid = uuidObj, mass = massObj,  size = sizeObj,  pos = posObj, pixel_size_um = resolutionObj, mean = meanSpcData, data = spectra)
  if(!is.null(posMotorsObj))
  {
    datacube$posMotors <- posMotorsObj
  }
  if(!is.null(normalizationsObj))
  {
    datacube$normalizations <- normalizationsObj
  }

  unlink(img_path, recursive = T)

  save(datacube, file = file.path(restore_path, "datacube.RImg"))
  if(!is.null(pb))
  {
    close(pb)
  }
  if(!is.null(fun_text))
  {
    fun_text("Done")
  }
  
  class(datacube) <- "rMSIObj"
  return(datacube)
}

#Re-Loads a previously loaded image which still have the ff files on HDD
#The LoadMsiData whill produce a R objecte image file in the same directory of ramdisk in HDD
#This function will check for the R objecte image in the provided path if it exist and ff object is correct it will return with the image object
#In case it is not possible loading the image it will return NULL
.FastLoad<-function( restorePath )
{
  #1- Check if restorePath exists
  if( !dir.exists(restorePath) )
  {
    cat("\nNo ramdisk directory, it will be created\n")
    return(NULL)
  }

  #2- Check if fast loading R image object file exists
  if( !file.exists( file.path(restorePath, "datacube.RImg") ))
  {
    cat("\nNo datacube.RImg in ramdisk, new ramdisk will be created\n")
    return(NULL)
  }

  #3- Load the Image object
  load(file.path(restorePath, "datacube.RImg") )

  #4- Check the ramdisk and load it
  ramDiskExists <- is.element(T, unlist(lapply(datacube$data, function(x) { file.exists(attr(attr(x, "physical"), "filename"))  })))
  if(!ramDiskExists)
  {
    cat("\nCurrent ramdisk has been corrupted, it will be created\n")
    return(NULL)
  }
  lapply(datacube$data, function(x){open(x)})

  cat("\nRamdisk has been sucessfully restored\n")
  return(datacube)
}


#'Remove an rMSI object ramdisk
#'
#' @param img an rMSI object.
#'
#' Removes a rMSI objecte ramdisk.
#'
#' @export
DeleteRamdisk<-function(img)
{
  lapply(img$data, function(x){ close(x) })
  ramdisk_path <- dirname(attr(attributes(img$data[[1]])$physical, "filename"))
  ramdisk_path_splited <- unlist(strsplit(ramdisk_path, "/"))
  if(.Platform$OS.type == "unix")
  {
    ramdisk_path <- "/"
  }
  else
  {
    ramdisk_path <- ""
  }
  ramdisk_path <- paste0(ramdisk_path,ramdisk_path_splited[1])
  if(length(ramdisk_path_splited) > 1)
  {
    for( i in 2:length(ramdisk_path_splited))
    {
      ramdisk_path<-file.path(ramdisk_path, ramdisk_path_splited[i])
      if( length(grep("ramdisk_", ramdisk_path_splited[i])) > 0 )
      {
        break;
      }
    }
  }
  unlink(ramdisk_path, recursive = T)
}

#' Create an empty rMSI object with defined mass axis and size.
#'
#' @param x_size the number of pixel in X direction.
#' @param y_size the number of pixel in Y direction.
#' @param mass_axis the mass axis.
#' @param pixel_resolution defined pixel size in um.
#' @param img_name the name for the image.
#' @param ramdisk_folder where ramdisk will be stored.
#' @param data_type a string determining data type used to store data.
#' @param uuid a string containing an universal unique identifier for the image. If it is not provided it will be created using a time code.
#'
#' Creates an empty rMSI object with the provided parameters. This method is usefull to implement importation of new data formats
#' and synthetic datasets to test and develop processing methods and tools.
#'
#' data_type possible values are:
#'  byte	  8 bit signed integer with NA.
#'  ubyte	  8 bit unsigned integer without NA.
#'  short   16 bit signed integer with NA.
#'  ushort  16 bit unsigned integer without NA.
#'  integer 32 bit signed integer with NA.
#'  single  32 bit float.
#'  double  64 bit float.
#'
#' @return the created rMSI object
#' @export
#'
CreateEmptyImage<-function(x_size, y_size, mass_axis, pixel_resolution, img_name = "New empty image", ramdisk_folder = getwd(), data_type = "integer", uuid = NULL)
{
  img<-list()
  img$name <- img_name
  img$mass <- mass_axis
  img$size <- c( x_size, y_size )
  names(img$size) <- c("x", "y")
  
  #Fill the UUID string
  if(is.null(uuid))
  {
    img$uuid <- uuid_timebased()
  }
  else
  {
    img$uuid <- uuid  
  }

  #Prepare the pos matrix
  img$pos <- matrix( ncol = 2, nrow = x_size*y_size )
  colnames(img$pos)<- c("x", "y")
  i <- 1
  for( xi in 1:x_size)
  {
    for( yi in 1:y_size)
    {
      img$pos[i,]<- c(xi, yi)
      i<-i+1
    }
  }

  img$pixel_size_um <-  pixel_resolution
  img$mean <- rep(0, length(mass_axis))

  #Prepare an empty datacube
  img$data<-.CreateEmptyRamdisk(length(mass_axis), nrow(img$pos), ramdisk_folder, vmode_type = data_type)

  class(img) <- "rMSIObj"
  return(img)
}

#' Create an empty rMSI object with defined mass axis and total number of pixels.
#'
#' @param num_of_pixels Total number of spectrums/pixels.
#' @param mass_axis the mass axis.
#' @param pixel_resolution defined pixel size in um.
#' @param img_name the name for the image.
#' @param ramdisk_folder where ramdisk will be stored.
#' @param data_type a string determining data type used to store data.
#' @param uuid a string containing an universal unique identifier for the image. If it is not provided it will be created using a time code.
#'
#' Creates an empty rMSI object with the provided parameters. This method is usefull to implement importation of new data formats
#' and synthetic datasets to test and develop processing methods and tools.
#' img$size is initialized with c(NA, NA) and the pos matrix with NA coords. Size and pos matrix must be filled by user.
#'
#' data_type possible values are:
#'  byte	  8 bit signed integer with NA.
#'  ubyte	  8 bit unsigned integer without NA.
#'  short   16 bit signed integer with NA.
#'  ushort  16 bit unsigned integer without NA.
#'  integer 32 bit signed integer with NA.
#'  single  32 bit float.
#'  double  64 bit float.
#'
#' @return the created rMSI object
#' @export
#'
CreateEmptyImage<-function(num_of_pixels, mass_axis, pixel_resolution, img_name = "New empty image", ramdisk_folder = getwd(), data_type = "integer", uuid = NULL)
{
  img<-list()
  img$name <- img_name
  img$mass <- mass_axis
  img$size <- c( NA, NA )
  names(img$size) <- c("x", "y")

  #Fill the UUID string
  if(is.null(uuid))
  {
    img$uuid <- format(Sys.time(), "%Y%m%d%H%M%S")
  }
  else
  {
    img$uuid <- uuid  
  }
  
  #Prepare the pos matrix
  img$pos <- matrix( NA, ncol = 2, nrow = num_of_pixels )
  colnames(img$pos)<- c("x", "y")

  img$pixel_size_um <-  pixel_resolution
  img$mean <- rep(0, length(mass_axis))

  #Prepare an empty datacube
  img$data<-.CreateEmptyRamdisk(length(mass_axis), nrow(img$pos), ramdisk_folder, vmode_type = data_type)
  class(img) <- "rMSIObj"
  return(img)
}

#' AverageSpectrum.
#' Computes the average spectrum of a whole rMSI image.
#'
#' @param img the rMSI image object.
#'
#' @return A vector with the average spectrum intensities. Masses are the same as the rMSI object.
#' @export
#'
  AverageSpectrum <- function(img)
{
  cat("Calculating Average Spectrum...\n")
  avgSpectrum <- tryCatch(
    {
      return(rMSIproc::AverageSpectrum_rMSIproc(img, NumOfThreads = min(parallel::detectCores()/2, 6)))
    },
    warning = function(war)
    {
        print(paste("WARNING in rMSI AverageSpectrum calling rMSIproc::AverageSpectrum_rMSIproc: ",war))
        return(NULL)
    }, 
    error = function(err) 
    {
      #No rMSIproc::AverageSpectrum_rMSIproc present... so using rMSI slow average
      pbavg <- txtProgressBar(min = 0, max = length(img$data), style = 3)
      avgI <- rep(0, length(img$mass))
      for( i in 1:length(img$data))
      {
        setTxtProgressBar(pbavg, i)
        avgI <- avgI +  colSums(img$data[[i]][,])
      }
      avgI <- avgI/nrow(img$pos)
      close(pbavg)
      return(avgI)
    })
  return(avgSpectrum)
}

#' BaseSpectrum.
#' Computes the base spectrum of a whole rMSI image where the intensity value
#' for each mass bin is calculated as the maxium of all mass channels.
#'
#' @param img the rMSI image object.
#'
#' @return A vector with the base spectrum intensities. Masses are the same as the rMSI object.
#' @export
#'
BaseSpectrum <- function(img)
{
  cat("Calculating Base Spectrum...\n")
  pb <- txtProgressBar(min = 0, max = length(img$data), style = 3)
  maxSub <- rep(0, length(img$mass))
  for( i in 1:length(img$data))
  {
    setTxtProgressBar(pb, i)
    maxSub <- maxSub + apply(img$data[[i]][,], 2, max)
  }
  close(pb)
  return(maxSub)
}

#' SortIDsByAcquisition: Order the rMSI pixel IDs according the acqusition sequence (first acquired pixel is the first one).
#'
#' @param img a rMSI image.
#'
#' @return a vector of ID's ordered acording acquisiton.
#' @export
#'
SortIDsByAcquisition <- function(img)
{
  idxArray <- matrix( c(1:nrow(img$pos), img$pos[,"x"], img$pos[,"y"]), nrow = nrow(img$pos), ncol = 3, byrow = F )
  colnames(idxArray) <- c("id", "x", "y")
  idxArray <- idxArray[order(idxArray[,"y"], idxArray[,"x"]), ]
  return(idxArray[,"id"])
}


#' PlotClusterImage.
#' Plot a segmentation image with the user-given clusters.
#'
#' @param posMat a two columns matrix where first column ara the x coodrinates of values and second column the y coordinates.
#' @param clusters a vector with integer number according the cluster of each pixel.
#' @param rotate rotation to apply.
#' @param pixel_size_um the pixel resolution in um.
#' @param labels_x x coordinates of text labels optionally overlaid to the plot.
#' @param labels_y y coordinates of text labels optionally overlaid to the plot.
#' @param labels_text text labels optionally overlaid to the plot.
#'
#' @return a vector with the used color for each cluster sorted according clustering numering in assending order.
#' @export
#'
PlotClusterImage <- function( posMat, clusters,  rotate = 0,  pixel_size_um = 100,
                              labels_x = NULL, labels_y = NULL, labels_text = NULL)
{
  img <- list()
  img$pos <- posMat
  colnames(img$pos) <- c("x", "y")
  img$size <- c(max(img$pos[ ,"x"]), max(img$pos[ ,"y"]))
  names(img$size) <- c("x", "y")
  img$pixel_size_um <- pixel_size_um

  #Prepare image matrix
  zplots<-matrix(0, nrow=img$size["x"], ncol=img$size["y"]) #Now I'm using a zero instead of NA to display a completely black background
  for( i in 1:nrow(img$pos))
  {
    zplots[img$pos[ i , "x" ], img$pos[ i , "y" ]] <- clusters[i]
  }

  #Create the raster
  my_raster <- raster::raster( nrow = ncol(zplots), ncol = nrow(zplots), xmn= 0, xmx= nrow(zplots), ymn= 0, ymx= ncol(zplots))
  raster::values(my_raster) <- as.vector(zplots)
  rm(zplots)

  #Put zplots matrix and some metadata in a list
  img_sgn <- list(raster = my_raster, mass = "", tolerance = 0, cal_resolution = img$pixel_size_um)
  rm(my_raster)

  oldPar <- par(no.readonly = T)
  raster_RGB<-.BuildSingleIonRGBImage(img_sgn, XResLevel = 3, light = 5)
  .plotMassImageRGB(raster_RGB, cal_um2pixels = img$pixel_size_um, rotation = rotate, display_axes = F, display_syscoords = F)

  #Get the colors used for clusters as a plotable form (RGB code) using the same rMSI internal function as raster image
  colras <- raster::raster( nrow = 1, ncol = 1+length(unique(clusters)))
  raster::values(colras) <- sort(c(0, unique(clusters)))
  rgbColRas <- .ReMappingIntensity2HSV(colras, value_multiplier = 5)
  clusterColors <- c()
  Rchannel <- raster::values(rgbColRas$layer.1)
  Gchannel <- raster::values(rgbColRas$layer.2)
  Bchannel <- raster::values(rgbColRas$layer.3)
  for( i in 1:length(unique(clusters))) #I'm avoiding the fist values because is the zero used to draw the background
  {
    clusterColors <- c(clusterColors, rgb( Rchannel[i + 1], Gchannel[i + 1], Bchannel[i + 1], 255, maxColorValue = 255))
  }

  text( x = labels_x, y = labels_y, labels = labels_text , adj = c(0.5, 0.0))
  par(oldPar)
  
  return(clusterColors)
}


#' PlotValues.
#'
#' Plot values in a image using the same methods as plotting an MS image.
#' The raster position of each value is definied in posMat.
#'
#' @param posMat a two columns matrix where first column ara the x coodrinates of values and second column the y coordinates.
#' @param values the values to plot.
#' @param rotate rotation to apply.
#' @param scale_title a text label for the color scale.
#' @param pixel_size_um the pixel resolution in um.
#' @param vlight the lighting of the plotted image.
#' @param labels_x x coordinates of text labels optionally overlaid to the plot.
#' @param labels_y y coordinates of text labels optionally overlaid to the plot.
#' @param labels_text text labels optionally overlaid to the plot.
#'
#' @export
#'
PlotValues <- function(posMat, values, rotate = 0, scale_title = "", pixel_size_um = 100, vlight = 5,
                       labels_x = NULL, labels_y = NULL, labels_text = NULL)
{
  fooImg <- list()
  fooImg$pos <- posMat
  colnames(fooImg$pos) <- c("x", "y")
  fooImg$size <- c(max(fooImg$pos[ ,"x"]), max(fooImg$pos[ ,"y"]))
  names(fooImg$size) <- c("x", "y")
  fooImg$pixel_size_um <- pixel_size_um
  PlotTICImage( fooImg, values, rotate, scale_title, vlight, labels_x = labels_x, labels_y = labels_y, labels_text = labels_text )
}

#' PlotTICImage.
#'
#' @param img an rMSI object.
#' @param TICs a vector of TIC values ordered acording pos array in img object.
#' @param rotate image rotation in degrees, possible values are: 0, 90, 180 and 270.
#' @param scale_title the title to show in scale axis (TIC by default).
#' @param vlight the lighting of the plotted image.
#' @param labels_x x coordinates of text labels optionally overlaid to the plot.
#' @param labels_y y coordinates of text labels optionally overlaid to the plot.
#' @param labels_text text labels optionally overlaid to the plot.
#'
#' @export
#'
PlotTICImage <- function(img, TICs = NULL, rotate = 0, scale_title = "TIC", vlight = 5, 
                         labels_x = NULL, labels_y = NULL, labels_text = NULL)
{
  #Calculate TICs
  if(is.null(TICs))
  {
    pb<-txtProgressBar(min = 0, max = length(img$data), style = 3 )
    setTxtProgressBar(pb, 0)
    TICs <- c()
    for( i in 1:length(img$data))
    {
      TICs <- c(TICs, rowSums(img$data[[i]][,]))
      setTxtProgressBar(pb, i)
    }
    close(pb)
  }

  #Do not plot infinites
  TICs[which(is.infinite(TICs))] <- 0

  #Prepare image matrix
  zplots<-matrix(0, nrow=img$size["x"], ncol=img$size["y"]) #Now I'm using a zero instead of NA to display a completely black background
  for( i in 1:nrow(img$pos))
  {
    zplots[img$pos[ i , "x" ], img$pos[ i , "y" ]] <- TICs[i]
  }

  #Create the raster
  my_raster <- raster::raster( nrow = ncol(zplots), ncol = nrow(zplots), xmn= 0, xmx= nrow(zplots), ymn= 0, ymx= ncol(zplots))
  raster::values(my_raster) <- as.vector(zplots)
  rm(zplots)

  #Put zplots matrix and some metadata in a list
  img_sgn <- list(raster = my_raster, mass = scale_title, tolerance = 0, cal_resolution = img$pixel_size_um)
  rm(my_raster)

  raster_RGB<-.BuildSingleIonRGBImage(img_sgn, XResLevel = 3, light = vlight)

  oldPar<- par(no.readonly = T)
  
  layout( matrix( 2:1, ncol = 2, nrow = 1, byrow = TRUE ), widths = c(7, rep(1, 2)) )
  .plotIntensityScale(img_sgn, light =  5)
  .plotMassImageRGB(raster_RGB, cal_um2pixels = img$pixel_size_um, rotation = rotate, display_axes = F, display_syscoords = F)
  
  text( x = labels_x, y = labels_y, labels = labels_text , adj = c(0.5, 0.0))
  
  par(oldPar)
}

#' ConvertrMSIimg2Bin.
#'
#' @param img a rMSI image object to be converted.
#' @param out_data_dir the resulting output data directory.
#'
#' @return nothing.
#' @export
#'
ConvertrMSIimg2Bin <- function( img, out_data_dir)
{

  dir.create(out_data_dir)

  #Export global m/z axis to txt
  write(img$mass, file.path(out_data_dir, "mz.txt" ), ncolumns = 1)

  #Export metadata as a plain ascii file:
  img_name_char <- paste("Original image name:\t", img$name, "\n", sep = "")
  img_size_charX <- paste("Image size in X:\t", img$size["x"], "\n", sep = "")
  img_size_charY <- paste("Image size in Y:\t", img$size["y"], "\n", sep = "")
  pixel_size_char <- paste("Pixel size (um):\t", img$pixel_size_um, "\n", sep = "")
  write( paste(img_name_char, img_size_charX, img_size_charY, pixel_size_char, sep =""), file= file.path(out_data_dir, "metadata.txt")  )

  #Export spectra intensity data
  spc_id <- 1
  subdir <- 1
  dir.create(file.path(out_data_dir, subdir))
  spc_in_subdir_count <- 0
  for( ic in 1:length(img$data))
  {
    cat(paste("Processing cube", ic, "of", length(img$data), "\n"))
    cube <- loadImgChunkFromCube(img, ic)

    for( i in 1:nrow(cube))
    {
      if(spc_in_subdir_count >= 1e4)
      {
        subdir <- subdir + 1
        dir.create(file.path(out_data_dir, subdir))
        spc_in_subdir_count <- 0
      }

      fname <- paste ("ID",spc_id, "_X", img$pos[spc_id, "x"] , "_Y", img$pos[spc_id, "y"], ".bin", sep = "")
      spc_id <- spc_id + 1
      spc <- cube[i, ]
      writeBin(spc, file.path(out_data_dir, subdir, fname ), size = 4, endian="little")
      spc_in_subdir_count <- spc_in_subdir_count + 1
    }

  }
}

#' ConvertBin2rMSIimg.
#'
#' @param in_data_dir data dir where the bin image is located.
#' @param out_img_tar_file if not NULL the imported image will be also stored as a tar file.
#'
#' @return the rMSI image object.
#' @export
#'
ConvertBin2rMSIimg <- function( in_data_dir, out_img_tar_file = NULL )
{
  #Get the metadata
  mz_axis <- as.numeric(read.table(file.path(in_data_dir, "mz.txt"))[,1])
  metadata <- read.table(file.path(in_data_dir, "metadata.txt"), sep = "\t", colClasses = "character")
  imgName <- as.character(metadata[1,2])
  xSize <- as.numeric(metadata[2,2])
  ySize <- as.numeric(metadata[3,2])
  resolution <- as.numeric(metadata[4,2])

  #List bin files
  bin_files<-list.files(path=in_data_dir, include.dirs = F, recursive = T, pattern = "*.bin", full.names = T)

  #Fill image fields
  img <-CreateEmptyImage( num_of_pixels = length(bin_files) , pixel_resolution = resolution, img_name =  imgName, mass_axis = mz_axis )
  img$size["x"] <- xSize
  img$size["y"] <- ySize

  #Prepare a dataframe with each bin file info
  print("Parsing bin file names...")
  pb<-txtProgressBar(min = 0, max = length(bin_files), style = 3 )
  pb_i <- 0
  ID <- c()

  for( bin_file in bin_files)
  {
    pb_i <- pb_i + 1
    setTxtProgressBar(pb, pb_i)
    pixel_fields <- strsplit(as.character(strsplit(basename(bin_file), split = "\\.")[[1]])[1], split = "_")[[1]]
    id <- as.numeric(strsplit(pixel_fields[1], split = "ID")[[1]])[2]
    X_cord <- as.numeric(strsplit(pixel_fields[2], split = "X")[[1]])[2]
    Y_cord <- as.numeric(strsplit(pixel_fields[3], split = "Y")[[1]])[2]
    img$pos[id, "x"] <- X_cord
    img$pos[id, "y"] <- Y_cord
    ID <- c(ID, id)
  }
  data_inf <- data.frame( ID, bin_files )
  data_inf <- data_inf[order(ID), ] #Sort by ID's (faster datacubes writing)
  close(pb)

  #Extract bin files
  print("Extracting spectra from bin files...")
  pb<-txtProgressBar(min = 0, max = nrow(data_inf), style = 3 )
  dm <- matrix(nrow = 100, ncol = length( mz_axis )) #Read HDD in chunks of 100 spectra
  partial_ids <- c()
  dm_irow <- 1
  for( i in 1:nrow(data_inf))
  {
    setTxtProgressBar(pb, i)
    dm[dm_irow, ] <- as.integer(readBin(as.character(data_inf[i, "bin_files"]), integer(),n=length(mz_axis) ,size = 4, signed=T, endian="little" ))
    partial_ids <- c(partial_ids, data_inf[i, "ID"])
    dm_irow <- dm_irow + 1

    if( dm_irow > nrow(dm) || i == nrow(data_inf))
    {
      saveImgChunkAtIds( img,  partial_ids, dm[1:length(partial_ids), ] )
      partial_ids <- c()
      dm_irow <- 1
    }
  }
  close(pb)

  print("Calculating average spectrum...")
  img$mean <- AverageSpectrum(img)

  if(!is.null(out_img_tar_file))
  {
    SaveMsiData(img, out_img_tar_file)
  }
  return(img)
}

#' plotParamAcqOrdered.
#'
#' @param img rMSI object from wich data must be ploted.
#' @param Param a vector of elements to plot.
#' @param yAxisLabel.
#'
#' Param will be ploted ordered according the order of pixels in MALDI acqusition.
#'
#' @export
#'
plotParamAcqOrdered <- function( img, Param, yAxisLabel = "Param" )
{
  idxArray <- matrix( c(1:nrow(img$pos), img$pos[,"x"], img$pos[,"y"]), nrow = nrow(img$pos), ncol = 3, byrow = F )
  colnames(idxArray) <- c("id", "x", "y")
  idxArray <- idxArray[order(idxArray[,"y"], idxArray[,"x"]), ]
  plot(Param[idxArray[,"id"]], type="l", col ="red", ylab = yAxisLabel, xlab = "Pixel" )
}



#' remap2ImageCoords.
#'
#' @param dataPos a pos matrix as it is in rMSI data object.
#'
#' This function should be only used to implement data importers from foreign formats.
#' This functions maps a MALDI motors coors space to a image coord space.
#' dataPos matrix is a two columns matrix where first column stores x positions and second y pixel positions.
#' a remapped dataPos matrix do not contain empty raster positions neighter offsets.
#'
#' @return the dataPos matrix remapped.
#' @export
#'
remap2ImageCoords <- function(dataPos)
{
  #1- Calc offsets and subtract it
  x_offset<-min(dataPos[,"x"])
  y_offset<-min(dataPos[,"y"])
  for(i in 1:nrow(dataPos))
  {
    dataPos[i, "x"] <- dataPos[i, "x"] - x_offset + 1
    dataPos[i, "y"] <- dataPos[i, "y"] - y_offset + 1
  }

  #2- Compute Motor coords range
  x_size<-max(dataPos[,"x"])
  y_size<-max(dataPos[,"y"])

  #3- Map MALDI motor coords to image cords (1-pixels steps)
  #It is important to map MALDI motors coords to image coords.
  #Otherwise, null extra pixels may be added leading to bad reconstruction
  px_map <- matrix( 0, nrow = x_size, ncol = y_size)
  for(i in 1:nrow(dataPos))
  {
    xi <- dataPos[i, "x"]
    yi <- dataPos[i, "y"]
    px_map[xi, yi]<- i
  }

  colNull <- which( base::colSums(px_map) == 0)
  rowNull <- which( base::rowSums(px_map) == 0)
  remap<-FALSE
  if( length(colNull) > 0 && length(rowNull) > 0 )
  {
    px_map_ <- px_map[ -rowNull , -colNull ]
    remap<-TRUE
  }
  if( length(colNull) > 0 && length(rowNull) == 0 )
  {
    px_map_ <- px_map[ , -colNull ]
    remap<-TRUE
  }
  if( length(colNull) == 0 && length(rowNull) > 0 )
  {
    px_map_ <- px_map[ -rowNull , ]
    remap<-TRUE
  }

  if(remap)
  {
    for(ix in 1:nrow(px_map_))
    {
      for(iy in 1:ncol(px_map_))
      {
        if(px_map_[ix, iy] > 0)
        {
          dataPos[px_map_[ix, iy], "x"] <- ix
          dataPos[px_map_[ix, iy], "y"] <- iy
        }
      }
    }
  }

  return(dataPos)
}

#' .controlled_loadAbort.
#'
#' @param text to be promt at R console
#' @param close_signal the function to call
#'
.controlled_loadAbort <- function(text, close_signal = NULL)
{
  if(!is.null(close_signal))
  {
    close_signal()
  }
  stop(text)
}

#' ParseBrukerXML.
#'
#' Reads a Bruker's xml file exported using fleximaging.
#' A list is returned where each element in the list is named according the ROI name.
#' Each element in the list consists in a data.frame with the pixels XY coordinates inside each ROI.
#'
#' @param xml_path the full path where XML file is stored.
#'
#' @return ROI pixel coordinates arranged in a named list.
#' @export
#' 
ParseBrukerXML <- function(xml_path)
{
  roilst <- CparseBrukerXML(path.expand(xml_path))
  
  if( !is.null(roilst$Error))
  {
    stop(roilst$Error)
  }
  
  return(roilst)
}

#' ReadBrukerRoiXML.
#' 
#' Reads a Bruker ROI XML file and matches it to an rMSI image object.
#' A list of rMSI ID's contained in each Bruker's ROI will be returned.
#'
#' @param img an rMSI MS image object.
#' @param xml_file a full path to a Bruker XML ROI file.
#'
#' @return a list containing lists of ID's for each ROI.
#' @export
#'
ReadBrukerRoiXML <- function(img, xml_file)
{
  if( is.null(img$posMotors))
  {
    stop("ERROR: image posMotros matrix not available.\nYou need to re-import MS data using a recent verison of rMSI.\n")
  }
  
  spectraListRois <- ParseBrukerXML(xml_file)
  lstRois <- list()
  imPosMat <- complex( real = img$posMotors[, "x"], imaginary = img$posMotors[, "y"])
  
  for( i in 1:length(spectraListRois))
  {
    cat(paste0("Parsing ROI ", i, " of ", length(spectraListRois), "\n"))
    lstRoisAux <- list(name = spectraListRois[[i]]$name, id = c())
    for( j in 1:nrow(spectraListRois[[i]]$pos))
    {
      #Extract original X Y Bruker Coords
      imCoord <- complex(real = spectraListRois[[i]]$pos$x[j], imaginary = spectraListRois[[i]]$pos$y[j])
      
      #Look for this Bruker coords in image pos matrix
      matchXY <- which(imPosMat == imCoord)
      if( length(matchXY) > 0)
      {
        lstRoisAux$id <- c( lstRoisAux$id, matchXY[1])
        if(  length(matchXY) > 1 )
        {
          cat(paste0("WARNING: roi ", spectraListRois[[i]]$name, " coordinates x", Re(imCoord), " , y", Im(imCoord), " are duplicated.\n" ))
        }
      }
    }
    
    if( length( lstRoisAux$id ) > 0 )
    {
      lstRoisAux$id <- sort( lstRoisAux$id, decreasing = F) #Sort ID's by assending order to get a faster disk access time
      lstRois[[(length(lstRois) + 1)]] <- lstRoisAux
    }
    else
    {
      cat(paste0("WARNING: deleting roi ",  spectraListRois[[i]]$name, " because it does not contain any pixel in the image.\n" ))
    }
  }
  
  return(lstRois)
}


#' CreateSubDataset.
#' 
#' Creates a new rMSI image object from sub-set of selected pixels by ID's.
#'
#' @param img the original rMSI object.
#' @param id a vector of ID's to retain in the sub data set.
#' @param ramdisk_path a full disk path where the new ramdisk will be stored.
#' @param new_mass a new mass axis if resampling must be used.
#'
#' @return the sub rMSI object.
#' @export
#'
CreateSubDataset <- function(img, id, ramdisk_path, new_mass = img$mass)
{
  cat("Creating the sub image...\n")

  if(!dir.exists(ramdisk_path))
  {
    dir.create(ramdisk_path, showWarnings = F, recursive = T)
  }
  
  #Resample data only if the supplied mass axis is different
  bResampleData <- !identical(img$mass, new_mass)
  
  subImg <- CreateEmptyImage(num_of_pixels = length(id), 
                             mass_axis = new_mass, 
                             pixel_resolution = img$pixel_size_um, 
                             img_name = paste0(img$name, "_sub"), 
                             ramdisk_folder = ramdisk_path, 
                             data_type = attr(attr(img$data[[1]], "physical"), "vmode"), 
                             uuid = img$uuid )
  
  subImg$pos <- remap2ImageCoords( img$pos[id, ] )
  
  if(!is.null(img$posMotors))
  {
    subImg$posMotors <- img$posMotors[id, ]
  }
  if( !is.null(img$normalizations))
  {
    subImg$normalizations <- img$normalizations
    for( i in 1:length(subImg$normalizations))
    {
      subImg$normalizations[[i]] <- subImg$normalizations[[i]][id]
    }
  }
  
  #Copy data on given id list
  pb <- txtProgressBar(min = 0, max = length(subImg$data), style = 3)
  istart <- 1
  for( i in 1:length(subImg$data))
  {
    setTxtProgressBar(pb, i)
    subID <- id[ istart:(istart + nrow(subImg$data[[i]]) - 1) ]
    istart <- istart +  nrow(subImg$data[[i]])
    dm <- loadImgChunkFromIds(img, subID)
    
    #Resampling...
    if(bResampleData)
    {
      dmSub <- matrix(nrow = nrow(dm), ncol = length(subImg$mass))
      for( irow in 1:nrow(dm))
      {
        dmSub[irow, ] <- approx(img$mass, dm[irow,], xout = subImg$mass, ties = "ordered", yleft = 0, yright = 0)$y
      }
    }
    else
    {
      dmSub <- dm
    }
    
    subImg$data[[i]][ , ] <- dmSub
  }
  close(pb)
  
  subImg$size <- c( max( subImg$pos[, "x"] ), max( subImg$pos[, "y"] ))
  names(subImg$size) <- c("x", "y")
  
  subImg$mean <- AverageSpectrum(subImg)
  return(subImg)
}

#' ROIAverageSpectra.
#' 
#' Calculates the average spectrum within each roi.
#'
#' @param img an rMSI object.
#' @param roi_list a roi list in the format returned by ReadBrukerRoiXML().
#'
#' @return all rois average spectra arranges in a list.
#' @export
#'
ROIAverageSpectra <- function( img, roi_list )
{
  if(length(roi_list) == 0)
  {
    stop("Error: emtpy region list\n")
  }
  
  cat("Calculating Average Spectra of ROI's...\n")

  roi_cubes <- lapply(roi_list, function(x){ getCubeRowFromIds(img, x$id) })
  rois_avg <- lapply(roi_list, function(x){list(name = x$name, mean = rep(0, length(img$mass)) )})
  
  pb <- txtProgressBar(min = 0, max = length(img$data), initial = 0, style = 3)
  LastId <- 0
  for( ic in 1:length(img$data) )
  {
    setTxtProgressBar(pb, ic)
    dm <- img$data[[ic]][,] #Load the complete data cube
    current_ids <- (LastId+1):(LastId+nrow(dm))
    
    for( ir in 1:length(roi_cubes))
    {
      idCube <- which(unlist(lapply(roi_cubes[[ir]], function(x){ x$cube})) == ic, arr.ind = T)
      if( length(idCube) > 0 )
      {
        idRows <- roi_cubes[[ir]][[idCube]]$row
        if(length(idRows) > 1)
        {
          rois_avg[[ir]]$mean <- rois_avg[[ir]]$mean + apply(dm[idRows,], 2, sum)
        }
        else
        {
          rois_avg[[ir]]$mean <- rois_avg[[ir]]$mean + dm[idRows,]
        }
      }
    }
    
    LastId <- current_ids[length(current_ids)]
  }
  
  #And divide by the pixel count on each roi
  for( ir in 1:length(rois_avg))
  {
    rois_avg[[ir]]$mean <- rois_avg[[ir]]$mean / length(roi_list[[ir]]$id)
  }
  close(pb)
  
  return(rois_avg)
}

#' ROIAverageSpectraByIds.
#' 
#' Calculates the average spectrum within each roi specified by a vector of pixel ID's
#'
#' @param img an rMSI object.
#' @param Ids Identifiers of spectra to use for average calculation.
#'
#' @return the ROI average spectrum.
#' @export
#'
ROIAverageSpectraByIds <- function( img, Ids )
{
  cat("Calculating Average Spectra of slected pixels...\n")
  
  roi_cubes <- getCubeRowFromIds(img, Ids)
  roi_avg <-  rep(0, length(img$mass))
  
  pb <- txtProgressBar(min = 0, max = length(img$data), initial = 0, style = 3)
  LastId <- 0
  for( ic in 1:length(img$data) )
  {
    setTxtProgressBar(pb, ic)
    dm <- img$data[[ic]][,] #Load the complete data cube
    current_ids <- (LastId+1):(LastId+nrow(dm))
    
    idCube <- which(unlist(lapply(roi_cubes, function(x){ x$cube})) == ic, arr.ind = T)
    if( length(idCube) > 0 )
    {
      idRows <- roi_cubes[[idCube]]$row
      if(length(idRows) > 1)
      {
        roi_avg <- roi_avg + apply(dm[idRows,], 2, sum)
      }
      else
      {
        roi_avg <- roi_avg + dm[idRows,]
      }
    }
    
    LastId <- current_ids[length(current_ids)]
  }
  
  #And divide by the pixel count
  roi_avg <- roi_avg / length(Ids)

  close(pb)
  
  return(roi_avg)
}

#' uuid.
#' 
#' Generates a timecode-based 16-bytes UUID.
#' The generated bytes are generated using the following pattern:
#' bytes 0 and 1: Year
#' byte 2: Month
#' byte 3: Day
#' byte 4: Hour
#' byte 5: Minute
#' byte 6: Second
#' bytes 7 to 15: random.
#'
#' @return the generated UUID encoded in a text string.
#'
uuid_timebased <- function()
{
  currentTime <- Sys.time()
  sUUID <- sprintf( "%04X", as.integer( format(currentTime, "%Y")))
  sUUID <- paste0( sUUID, sprintf( "%02X", as.integer( format(currentTime, "%m"))))
  sUUID <- paste0( sUUID, sprintf( "%02X", as.integer( format(currentTime, "%d"))))
  sUUID <- paste0( sUUID, sprintf( "%02X", as.integer( format(currentTime, "%H"))))
  sUUID <- paste0( sUUID, sprintf( "%02X", as.integer( format(currentTime, "%M"))))
  sUUID <- paste0( sUUID, sprintf( "%02X", as.integer( format(currentTime, "%S"))))
  sUUID <-paste0(sUUID, paste0(sprintf("%02X",sample(0:255, 9)), collapse = ""))
  
  return(sUUID)
}