R/display_3D_stack.R

In espadon: Easy Study of Patient DICOM Data in Oncology

#' Display in 3D the selected planes of an \pkg{espadon} class volume
#' @description The \code{display.3D.stack} function displays in 3D the requested 
#' cutting planes of a "volume" class object.
#' @param vol "volume" class object to display.
#' @param k.idx vector of cutting plane numbers to be displayed, to be chosen in 
#' \code{vol$k.idx}. By default \code{k.idx} is a vector of 10 uniformly 
#' distributed cutting planes in the volume.
#' @param display.ref Character string. Pseudonym of the frame of reference used 
#' for display.
#' @param T.MAT "t.mat" class object, created by \link[espadon]{load.patient.from.Rdcm} 
#' or \link[espadon]{load.T.MAT}. If \code{T.MAT} is \code{NULL}, \code{vol} must 
#' be displayed in \code{display.ref = vol$ref.pseudo}.
#' @param col Vector, representing the color palette of the image. Transparent 
#' colors are not represented. 
#' @param breaks  One of :
#' \itemize{
#' \item \code{NULL} : The minimum and the maximum value of the \code{vol} define 
#' the range.
#' \item Vector giving the breakpoints of each color.
#' }
#' @param cube Boolean. If \code{TRUE} (default), the "volume" edges are displayed.
#' @param border Boolean. If \code{TRUE} (default), the borders of the planes defined 
#' in \code{k.idx} are displayed.
#' @param ktext Boolean. If \code{TRUE} (default), the  selected cutting plane numbers 
#' are displayed.
#' @param line.col Color of cube, planes and texts displayed.
#' @param line.lwd Line width of the border and cube, by default at 1. 
#' @param cex Numeric character expansion factor of displayed plan numbers.

#' @return  Returns a display of the \code{k.idx} cutting planes of \code{vol}, 
#' in the current \pkg{RGL} window if it exists, in a new window otherwise. The 
#' colors of the palettes are managed by \code{col} and \code{breaks}.  

#' @examples
#' # loading of toy-patient objects (decrease dxyz for better result)
#' step <- 4
#' patient <- toy.load.patient (modality = "ct", dxyz = rep (step, 3))
#' 
#' # display o 3 planes
#' rgl::open3d()
#' display.3D.stack (patient$ct[[1]],
#'                   col = pal.RVV (200, alpha = c(rep(0,90), rep (1, 110))))



#' @export
#' @importFrom rgl mesh3d wire3d segments3d texts3d material3d
#' @importFrom methods is
display.3D.stack <- function (vol,  
                             k.idx = unique(vol$k.idx[seq(1,vol$n.ijk[3], length.out=10)]), 
                             display.ref = vol$ref.pseudo, 
                             T.MAT = NULL,  
                             col = grey.colors(10, start = 0, end = 1, alpha = c(rep(0,1),rep(1,9))),
                             breaks = NULL,
                             cube = TRUE, 
                             border = TRUE,
                             ktext = TRUE,
                             line.col = "#379DA2",
                             line.lwd = 1,
                             cex = 1){
  
  if (!is (vol, "volume")) 
    stop ("vol should be a volume class object.")
  
  if (is.null(vol$vol3D.data)) 
    stop  ("empty vol$vol3D.data.")
  
  M <- get.rigid.M (T.MAT, vol$ref.pseudo, display.ref)
  if (is.null (M)) stop ("cannot display anything in selected frame of reference.")
  
  imin <- as.numeric(vol$min.pixel)
  imax <- as.numeric(vol$max.pixel) 
  
  if (!is.null(breaks) & !any(is.na(breaks))) {
    if (length(breaks)!= length (col) + 1) stop ("length(breaks) must be equal to length(col) + 1.")
    
  } else {
    breaks <- .pixel.scale (imin,imax,length(col))
  }
  if (breaks[1]>imin) breaks[1] <- imin - 1
  if (breaks[length (breaks)]<imax) breaks[length (breaks)] <- imax + 1
  
  color  <- substr(col,1,7)
  alpha <- substr(col,8,9)
  if (nchar(alpha[1])==0) alpha= rep ("ff",length(col))
  cube_disp <- vol$cube.idx
  cube_disp[c(1:3,6:7,11,13,15,17:18,22,29)] <- cube_disp[c(1:3,6:7,11,13,15,17:18,22,29)] - 0.5
  cube_disp[c(5,9:10,14,19,21,23,25:27,30:31)] <- cube_disp[c(5,9:10,14,19,21,23,25:27,30:31)] + 0.5
  
  if (cube & !any (apply(vol$cube.idx == 0 ,1, all))){

    line <-unique (do.call(rbind.data.frame,lapply(1:ncol(vol$cube.idx), function(i) 
      data.frame(a=rep(i,3),b=which(apply(sweep(vol$cube.idx,1,vol$cube.idx[,i])==0,2,sum)==3)))))

    pt <- (M %*% vol$xyz.from.ijk %*% cube_disp)[1:3,]
    for (i in 1:nrow(line)) rgl::segments3d (t (pt[ ,c (line[i,1], line[i,2])]), color = line.col, lwd=line.lwd)
  }

  for (k in k.idx){
    p <- cube_disp[,1:4]
    p[3,] <- k
    pt <- (M %*% vol$xyz.from.ijk %*% p)[1:3,]
    if (border){    
      segments3d (t (pt[,c(1,2)]), color=line.col, lwd=line.lwd)
      segments3d (t (pt[,c(2,3)]), color=line.col, lwd=line.lwd)
      segments3d (t (pt[,c(3,4)]), color=line.col, lwd=line.lwd)
      segments3d (t (pt[,c(4,1)]), color=line.col, lwd=line.lwd)
      
      if (ktext) texts3d(t (pt[,c(1)]), texts=paste0("k=",k," "), 
                         color=line.col,cex=cex, adj=1)  
    }
      
    idx <- which (k == vol$k.idx)
    map <-  array(vol$vol3D.data[ , ,idx],dim=vol$n.ijk[1:2])
    
    if(!all(is.na(map))){
      cut.map <- cut (as.numeric (map),breaks, include.lowest=TRUE)
      alpha.layer <-  matrix (alpha[cut.map], nrow=nrow(map))
      map.layer <- matrix (color[cut.map], nrow=nrow(map))
      ind <- which(alpha.layer!="00", arr.ind = T) 
      
      if(nrow(ind)>0){
        pt.ijk <- as.matrix(cbind(as.data.frame(ind)-1,k=vol$k.idx[idx],t=1))
        
        pt <- (as.matrix(cbind(as.data.frame(ind)-1,k=vol$k.idx[idx],t=1)) %*% t(M %*% vol$xyz.from.ijk))
        A <- (as.matrix (sweep(pt.ijk,2, c ( 0.5,  0.5, 0, 0))) %*% t(M %*% vol$xyz.from.ijk))
        B <- (as.matrix (sweep(pt.ijk,2, c (-0.5,  0.5, 0, 0))) %*% t(M %*% vol$xyz.from.ijk))
        C <- (as.matrix (sweep(pt.ijk,2, c (-0.5, -0.5, 0, 0))) %*% t(M %*% vol$xyz.from.ijk))
        D <- (as.matrix (sweep(pt.ijk,2, c ( 0.5, -0.5, 0, 0))) %*% t(M %*% vol$xyz.from.ijk))
        
        mesh <- list()
        mesh$vb <- t((rbind(A,B,C,D)))
        mesh$qu <- rbind(1:nrow(A),(1:nrow(A))+nrow(A),(1:nrow(A))+2*nrow(A),(1:nrow(A))+3*nrow(A))
        material <-material3d()
        material$specular = "black"
        material$color <- as.character(map.layer[ind])
        material$lit <- FALSE
        m <- mesh3d(vertices=mesh$vb, quads=mesh$qu , meshColor="faces", material=material)
        
        shade3d(m)
      }
    }
  }

}