R/glrlm.R

In neuroconductor-devel-releases/RIA: Radiomics Image Analysis Toolbox for Medial Images

#' @title Creates gray-level run length matrix from RIA image
#' @export
#' @encoding UTF-8
#'
#' @description  Creates gray-level run length matrix (GLRLM) from \emph{RIA_image}.
#' GLRLM assesses the spatial relation of voxels to each other by investigating how many times
#' same value voxels occur next to each other in a given direction. By default the \emph{$modif}
#' image will be used to calculate GLRLMs. If \emph{use_slot} is given, then the data
#' present in \emph{RIA_image$use_slot} will be used for calculations.
#' Results will be saved into the \emph{glrlm} slot. The name of the subslot is determined
#' by the supplied string in \emph{save_name}, or is automatically generated by RIA. \emph{off_right},
#' \emph{off_down} and \emph{off_z} logicals are used to indicate the direction of the runs.
#'
#' @param RIA_data_in \emph{RIA_image}.
#'
#' @param off_right integer, positive values indicate to look to the right, negative values
#' indicate to look to the left, while 0 indicates no offset in the X plane.
#'
#' @param off_down integer, positive values indicate to look to the right, negative values
#' indicate to look to the left, while 0 indicates no offset in the Y plane.
#'
#' @param off_z integer, positive values indicate to look to the right, negative values
#' indicate to look to the left, while 0 indicates no offset in the Z plane.
#'
#' @param use_type string, can be \emph{"single"} which runs the function on a single image,
#' which is determined using \emph{"use_orig"} or \emph{"use_slot"}. \emph{"discretized"}
#' takes all datasets in the \emph{RIA_image$discretized} slot and runs the analysis on them.
#'
#' @param use_orig logical, indicating to use image present in \emph{RIA_data$orig}.
#' If FALSE, the modified image will be used stored in \emph{RIA_data$modif}.
#'
#' @param use_slot string, name of slot where data wished to be used is. Use if the desired image
#' is not in the \emph{data$orig} or \emph{data$modif} slot of the \emph{RIA_image}. For example,
#' if the desired dataset is in \emph{RIA_image$discretized$ep_4}, then \emph{use_slot} should be
#' \emph{discretized$ep_4}. The results are automatically saved. If the results are not saved to
#' the desired slot, then please use \emph{save_name} parameter.
#'
#' @param save_name string, indicating the name of subslot of \emph{$glcm} to save results to.
#' If left empty, then it will be automatically determined based on the
#' last entry of \emph{RIA_image$log$events}.
#'
#' @param verbose_in logical indicating whether to print detailed information.
#' Most prints can also be suppressed using the \code{\link{suppressMessages}} function.
#'
#' @return \emph{RIA_image} containing the GLRLM.
#'
#' @examples \dontrun{
#' #Discretize loaded image and then calculate GLRLM matrix of RIA_image$modif
#' RIA_image <- discretize(RIA_image, bins_in = c(4, 8), equal_prob = TRUE,
#' use_orig = TRUE, write_orig = FALSE)
#' RIA_image <- glrlm(RIA_image, use_orig = FALSE, verbose_in = TRUE)
#'
#' #Use use_slot parameter to set which image to use
#' RIA_image <- glrlm(RIA_image, use_orig = FALSE, use_slot = "discretized$ep_4",
#' off_right = 1, off_down = 1, off_z = 0)
#' 
#' #Batch calculation of GLRLM matrices on all discretized images
#' RIA_image <- glrlm(RIA_image, use_type = "discretized",
#' off_right = 1, off_down = 1, off_z = 0)
#' }
#' 
#' @references 
#' Mary M. Galloway et al.
#' Texture analysis using gray level run lengths.
#' Computer Graphics and Image Processing. 1975; 4:172-179.
#' DOI: 10.1016/S0146-664X(75)80008-6
#' \url{https://www.sciencedirect.com/science/article/pii/S0146664X75800086}
#' 
#' Márton KOLOSSVÁRY et al.
#' Radiomic Features Are Superior to Conventional Quantitative Computed Tomographic
#' Metrics to Identify Coronary Plaques With Napkin-Ring Sign
#' Circulation: Cardiovascular Imaging (2017).
#' DOI: 10.1161/circimaging.117.006843
#' \url{https://www.ncbi.nlm.nih.gov/pubmed/29233836}
#' 
#' Márton KOLOSSVÁRY et al.
#' Cardiac Computed Tomography Radiomics: A Comprehensive Review on Radiomic Techniques.
#' Journal of Thoracic Imaging (2018).
#' DOI: 10.1097/RTI.0000000000000268
#' \url{https://www.ncbi.nlm.nih.gov/pubmed/28346329}
#' @encoding UTF-8

glrlm <- function(RIA_data_in, off_right = 1, off_down = 0, off_z = 0, use_type = "single", use_orig = FALSE, use_slot = NULL, save_name = NULL, verbose_in = TRUE)
{
    data_in_orig <- check_data_in(RIA_data_in, use_type = use_type, use_orig = use_orig, verbose_in = verbose_in)
    
    if(any(class(data_in_orig) != "list")) data_in_orig <- list(data_in_orig)
    list_names <- names(data_in_orig)
    if(!is.null(save_name) & (length(data_in_orig) != length(save_name))) {stop(paste0("PLEASE PROVIDE THE SAME NUMBER OF NAMES AS THERE ARE IMAGES!\n",
                                                                                       "NUMBER OF NAMES:  ", length(save_name), "\n",
                                                                                       "NUMBER OF IMAGES: ", length(data_in), "\n"))
    }
    
    for (k in 1: length(data_in_orig))
    {
        data_in <-  data_in_orig[[k]]
        
        if(off_z & dim(data_in)[3] == 1) {{stop("WARNING: CANNOT ASSESS Z PLANE OFFSET IF DATA IS 2D!")}}
        
        data_NA <- as.vector(data_in)
        data_NA <- data_NA[!is.na(data_NA)]
        if(length(data_NA) == 0) {stop("WARNING: SUPPLIED RIA_image DOES NOT CONTAIN ANY DATA!!!")}
        if(length(dim(data_in)) < 2 | length(dim(data_in)) > 3) stop(paste0("DATA LOADED IS ", length(dim(data_in)), " DIMENSIONAL. ONLY 2D AND 3D DATA ARE SUPPORTED!"))
        
        
        dim_x <- dim(data_in)[1]
        dim_y <- dim(data_in)[2]
        dim_z <- ifelse(!is.na(dim(data_in)[3]), dim(data_in)[3], 1)
        
        if(off_right > 0) off_right = 1; if(off_right < 0) off_right = -1
        if(off_down > 0) off_down = 1; if(off_down < 0) off_down = -1
        if(off_z > 0) off_z = 1; if(off_z < 0) off_z = -1
        
        base_m <- array(NA, dim = c(dim_x+dim_x*abs(off_down), dim_y+dim_y*abs(off_right), dim_z+dim_z*abs(off_z)))
        offset <- array(c(dim_x*off_down, dim_y*off_right, dim_z*off_z)); offset[offset == 0] <- NA; offset <- min(abs(offset), na.rm = TRUE)
        
        #Position data into correct possition
        if(off_down > -1 & off_right > -1 & off_z > -1)          {base_m[1:dim_x, 1:dim_y, 1:dim_z] <- data_in
        }else if(off_down == -1 & off_right > -1 & off_z > -1)   {base_m[(dim_x+1):(2*dim_x), 1:dim_y, 1:dim_z] <- data_in
        }else if(off_down > -1 & off_right == -1 & off_z > -1)   {base_m[1:dim_x, (dim_y+1):(2*dim_y), 1:dim_z] <- data_in
        }else if(off_down > -1 & off_right > -1 & off_z == -1)   {base_m[1:dim_x, 1:dim_y, (dim_z+1):(2*dim_z)] <- data_in
        }else if(off_down == -1 & off_right == -1 & off_z > -1)  {base_m[(dim_x+1):(2*dim_x), (dim_y+1):(2*dim_y), 1:dim_z] <- data_in
        }else if(off_down == -1 & off_right > -1 & off_z == -1)  {base_m[(dim_x+1):(2*dim_x), 1:dim_y, (dim_z+1):(2*dim_z)] <- data_in
        }else if(off_down > -1 & off_right == -1 & off_z == -1)  {base_m[1:dim_x, (dim_y+1):(2*dim_y), (dim_z+1):(2*dim_z)] <- data_in
        }else if(off_down == -1 & off_right == -1 & off_z == -1) {base_m[(dim_x+1):(2*dim_x), (dim_y+1):(2*dim_y), (dim_z+1):(2*dim_z)] <- data_in
        }else {stop("WARNING: OFFSETS ARE NOT APPROPRIATE PLEASE SEE help(glrlm)!!!")
        }
        
        #create gray level number, first by the name of the file, then the event log, then by the number of gray levels
        num_ind <- unlist(gregexpr('[1-9]', list_names[k]))
        num_txt <- substr(list_names[k], num_ind[1], num_ind[length(num_ind)])
        gray_levels <- as.numeric(num_txt)
        if(length(gray_levels) == 0) {
            txt <- automatic_name(RIA_data_in, use_orig, use_slot)
            num_ind <- unlist(gregexpr('[1-9]', txt))
            num_txt <- substr(txt, num_ind[1], num_ind[length(num_ind)])
            gray_levels <- as.numeric(num_txt)
        } else {data_v <- as.vector(data_in); data_v <- data_v[!is.na(data_v)]; gray_levels <- length(unique(data_v))
        }
        
        glrlm <- array(NA, c(gray_levels, offset))
        
        for (i in 1: gray_levels) {
            
            base_filt_m <- data_in; base_filt_m[base_filt_m != i] <- NA; base_filt_m[base_filt_m == i] <- 1
            base_filt_change_m <- array(NA, dim = c(dim_x+dim_x*abs(off_down), dim_y+dim_y*abs(off_right), dim_z+dim_z*abs(off_z)))
            
            if(off_down > -1 & off_right > -1 & off_z > -1)           {base_filt_change_m[1:dim_x, 1:dim_y, 1:dim_z] <- base_filt_m
            } else if(off_down == -1 & off_right > -1 & off_z > -1)   {base_filt_change_m[(dim_x+1):(2*dim_x), 1:dim_y, 1:dim_z] <- base_filt_m
            } else if(off_down > -1 & off_right == -1 & off_z > -1)   {base_filt_change_m[1:dim_x, (dim_y+1):(2*dim_y), 1:dim_z] <- base_filt_m
            } else if(off_down > -1 & off_right > -1 & off_z == -1)   {base_filt_change_m[1:dim_x, 1:dim_y, (dim_z+1):(2*dim_z)] <- base_filt_m
            } else if(off_down == -1 & off_right == -1 & off_z > -1)  {base_filt_change_m[(dim_x+1):(2*dim_x), (dim_y+1):(2*dim_y), 1:dim_z] <- base_filt_m
            } else if(off_down == -1 & off_right > -1 & off_z == -1)  {base_filt_change_m[(dim_x+1):(2*dim_x), 1:dim_y, (dim_z+1):(2*dim_z)] <- base_filt_m
            } else if(off_down > -1 & off_right == -1 & off_z == -1)  {base_filt_change_m[1:dim_x, (dim_y+1):(2*dim_y), (dim_z+1):(2*dim_z)] <- base_filt_m
            } else if(off_down == -1 & off_right == -1 & off_z == -1) {base_filt_change_m[(dim_x+1):(2*dim_x), (dim_y+1):(2*dim_y), (dim_z+1):(2*dim_z)] <- base_filt_m
            } else {stop("WARNING: OFFSETS ARE NOT APPROPRIATE PLEASE SEE help(glrlm)!!!")
            }
            
            for (j in 1: (offset-1)) {
                shift_m <- array(NA, dim = c(dim_x+dim_x*abs(off_down), dim_y+dim_y*abs(off_right), dim_z+dim_z*abs(off_z)))
                
                if(off_down > -1 & off_right > -1 & off_z > -1)           {shift_m[(1+j*off_down):(dim_x+j*off_down), (1+j*off_right):(dim_y+j*off_right), (1+j*off_z):(dim_z+j*off_z)] <- base_filt_m
                } else if(off_down == -1 & off_right > -1 & off_z > -1)   {shift_m[((dim_x+1)+j*off_down):((2*dim_x)+j*off_down), (1+j*off_right):(dim_y+j*off_right), (1+j*off_z):(dim_z+j*off_z)] <- base_filt_m
                } else if(off_down > -1 & off_right == -1 & off_z > -1)   {shift_m[(1+j*off_down):(dim_x+j*off_down), ((dim_y+1)+j*off_right):((2*dim_y)+j*off_right), (1+j*off_z):(dim_z+j*off_z)] <- base_filt_m
                } else if(off_down > -1 & off_right > -1 & off_z == -1)   {shift_m[(1+j*off_down):(dim_x+j*off_down), (1+j*off_right):(dim_y+j*off_right), ((dim_z+1)+j*off_z):((2*dim_z)+j*off_z)] <- base_filt_m
                } else if(off_down == -1 & off_right == -1 & off_z > -1)  {shift_m[((dim_x+1)+j*off_down):((2*dim_x)+j*off_down), ((dim_y+1)+j*off_right):((2*dim_y)+j*off_right), (1+j*off_z):(dim_z+j*off_z)] <- base_filt_m
                } else if(off_down == -1 & off_right > -1 & off_z == -1)  {shift_m[((dim_x+1)+j*off_down):((2*dim_x)+j*off_down), (1+j*off_right):(dim_y+j*off_right), ((dim_z+1)+j*off_z):((2*dim_z)+j*off_z)] <- base_filt_m
                } else if(off_down > -1 & off_right == -1 & off_z == -1)  {shift_m[(1+j*off_down):(dim_x+j*off_down), ((dim_y+1)+j*off_right):((2*dim_y)+j*off_right), ((dim_z+1)+j*off_z):((2*dim_z)+j*off_z)] <- base_filt_m
                } else if(off_down == -1 & off_right == -1 & off_z == -1) {shift_m[((dim_x+1)+j*off_down):((2*dim_x)+j*off_down), ((dim_y+1)+j*off_right):((2*dim_y)+j*off_right), ((dim_z+1)+j*off_z):((2*dim_z)+j*off_z)] <- base_filt_m
                } else {stop("WARNING: OFFSETS ARE NOT APPROPRIATE PLEASE SEE help(glrlm)!!!")
                }
                
                diff_m <- base_filt_change_m - shift_m
                count_diff <- length(diff_m[!is.na(diff_m)])
                glrlm[i,(j+1)] <- count_diff
                
                base_filt_change_m <- diff_m
            }
            
            count_gl <- base_filt_m; count_gl <- count_gl[!is.na(count_gl)]; count_gl <- sum(count_gl, na.rm = TRUE)
            
            for (p in seq(dim(glrlm)[2], 3, -1)) {
                if(glrlm[i, p] > 0) {
                    
                    m = 2
                    for (q in seq((p-1), 2, -1)) {
                        glrlm[i, q] <- glrlm[i, q] - glrlm[i, p]*m
                        m = m+1
                    }
                }
            }
            
            glrlm_r_sum <- sum(col(glrlm)[1,]*glrlm[i,], na.rm = TRUE)
            if(glrlm_r_sum != count_gl) {glrlm[i,1] <- (count_gl-glrlm_r_sum)
            } else {glrlm[i,1] <- 0}
            
        }
        
        if(use_type == "single") {
            
            if(any(class(RIA_data_in) == "RIA_image") )
            {
                if(is.null(save_name)) {
                    txt <- automatic_name(RIA_data_in, use_orig, use_slot)
                    txt <- paste0(txt, "_", as.numeric(off_right), as.numeric(off_down), as.numeric(off_z))
                    
                    RIA_data_in$glrlm[[txt]] <- glrlm
                    
                }
                if(!is.null(save_name)) {RIA_data_in$glrlm[[save_name]] <- glrlm
                }
            }
        }
        
        if(use_type == "discretized") {
            if(any(class(RIA_data_in) == "RIA_image"))
            {
                if(is.null(save_name[k])) {
                    txt <- list_names[k]
                    txt <- paste0(txt, "_", as.numeric(off_right), as.numeric(off_down), as.numeric(off_z))
                    
                    RIA_data_in$glrlm[[txt]] <- glrlm
                }
                if(!is.null(save_name[k])) {RIA_data_in$glrlm[[save_name[k]]] <- glrlm
                }
            }
        }
        
        if(is.null(save_name)) {txt_name <- txt
        } else {txt_name <- save_name}
        if(verbose_in) {message(paste0("GLRLM WAS SUCCESSFULLY ADDED TO '", txt_name, "' SLOT OF RIA_image$glrlm\n"))}
        
        
    }
    
    if(any(class(RIA_data_in) == "RIA_image")) {return(RIA_data_in)
    } else {return(glrlm)}
    
}