R/toy_rtdose.R

In espadon: Easy Study of Patient DICOM Data in Oncology

#' @importFrom stats pnorm 
.toy.rtdose.from.roi <- function(D.max, DSA, beam.nb, vol, struct, roi.name = NULL, 
                                 roi.sname = NULL, roi.idx = NULL, alias = "", 
                                 description = "") {
  
  roi.idx <- select.names(struct$roi.info$roi.pseudo, roi.name, roi.sname, roi.idx)
  ct.ptv <- nesting.roi(vol, struct, roi.idx = roi.idx, 
                         xyz.margin = rep(4 * max(vol$dxyz),3))
  bin.PTV <- bin.from.roi(ct.ptv,struct, roi.sname = "ptv", verbose = FALSE)
  bin.PTV.d <- bin.dilation(bin.PTV, radius = 2 * max(abs(vol$dxyz)))

  PTV.G <- as.numeric(struct$roi.info[roi.idx, c("Gx","Gy","Gz")])
  alpha <- atan(max(c(abs(as.numeric(struct$roi.info[roi.idx, c("max.x","max.y","max.z")] - 
                                   PTV.G)) + abs(vol$dxyz),
                      abs(as.numeric(struct$roi.info[roi.idx, c("min.x","min.y","min.z")] - 
                                   PTV.G)) + abs(vol$dxyz)))/DSA)*180/pi

  theta <- rev(seq(360, 0, length.out = beam.nb + 1)[-1] * pi / 180)
  
  direction <- cbind(-sin(theta), cos(theta), 0)
  orientation <- as.matrix(do.call(rbind,
                                   lapply(1:beam.nb, function(i) 
                                      c(c(0,0,1), vector.product(direction[i,],c(0,0,1))))))
  src <- sweep(-DSA * direction, 2, PTV.G,FUN = "+")
  
  src.extrem.pt.dist <- max(apply(src,1, function(sc)  
    max(.fnorm(sweep(t((vol$xyz.from.ijk %*% vol$cube.idx)[1:3,]), 2, sc)))))
  dalpha <- 180 * min(vol$dxyz) / src.extrem.pt.dist / pi
  alpha <- alpha + dalpha
  
  density <- vol.copy(vol)
  density$vol3D.data[density$vol3D.dat <= -900]  <- -1000
  density$vol3D.data[]  <- (density$vol3D.data[] / 1000 + 1)
  density$max.pixel <- (density$max.pixel / 1000 + 1)
  density$min.pixel <- (density$min.pixel / 1000 + 1)
  
  D.list <- lapply(1:beam.nb, function(theta.idx){
    
    fan <- fan.beam(alpha = alpha,dalpha, origin = as.numeric(src[theta.idx,]),
                    orientation =  as.numeric(orientation[theta.idx,]), 
                    ref.pseudo = bin.PTV.d$ref.pseudo)
    vox <- fan.to.voxel(bin.PTV.d,fan,restrict = T )
    ray.index <- sort( unique(vox$ray.index))
    fan$xyz <- fan$xyz[ray.index, ]
    fan$local.coord <- fan$local.coord[ray.index, ]
    r <- as.numeric(DSA * sqrt((fan$xyz %*% fan$orientation[1:3])^2 + (fan$xyz %*% fan$orientation[4:6])^2))
    flou <- 1 - pnorm(r,tan(alpha*pi/180) * DSA * 0.9 ,1)
    voxel <- .fan.to.voxel.with.att(density,fan, att = TRUE)
    voxel$att  <-  pnorm(voxel$att, 2, 5) * exp(-voxel$att / 150)
    vol.idx <- unique(voxel$vol.index)
    D_ <- vol.copy(vol, modality = "rtdose")
    D_$vol3D.data[] <- 0
    D_$vol3D.data <-  .mean_voxC(as.numeric(D_$vol3D.data), vol.idx - 1, 
                                 voxel$vol.index - 1, voxel$att, flou[voxel$ray.index])
    D_$vol3D.data  <- array(D_$vol3D.data ,dim = vol$n.ijk)
    D_$vol3D.data  <- D_$vol3D.data / max(D_$vol3D.data,na.rm = TRUE)
    D_$max.pixel <- 1
    D_$min.pixel <- 0
    D_
  })
  
  D <- vol.copy(vol,modality = "rtdose" )
  D$vol3D.data[] <- 0
  D$min.pixel <- 0
  D$max.pixel <- 0
  for (i in 1:beam.nb) D <- vol.sum(D, D.list[[i]], alias = alias, description = description)
  
  D$vol3D.data <- D.max*D$vol3D.data / nesting.roi(D, struct, 
                                                   roi.idx = roi.idx)$max.pixel
  D$max.pixel  <- max(D$vol3D.data,na.rm = TRUE)
  D$min.pixel  <- min(D$vol3D.data,na.rm = TRUE)
  return(D)
}

.toy.rtdose.from.bin <- function(D.max, DSA, beam.nb, vol, bin , alias = "", description = "") {
  
  PTV.xyz <- get.xyz.from.index(which(bin$vol3D.data),bin)
  margin <- max(abs(vol$dxyz))
  pt.min <- apply(PTV.xyz,2,min)
  pt.max <- apply(PTV.xyz,2,max)
  
  bin.PTV.d <- nesting.cube(bin, pt.min = pt.min - 4 * margin,
                            pt.max =  pt.max + 4 * margin)
  
  PTV.G <-  apply(PTV.xyz,2,mean)
  alpha <- atan(max(c(abs(as.numeric(pt.max - PTV.G)) + abs(vol$dxyz),
                      abs(as.numeric(pt.min - PTV.G)) + abs(vol$dxyz))) / DSA) * 180 / pi
  theta <- rev(seq(360, 0, length.out = beam.nb + 1)[-1] * pi / 180)
  
  direction <- cbind(-sin(theta), cos(theta), 0)
  orientation <- as.matrix(do.call(rbind,
                                   lapply(1:beam.nb, function(i) 
                                      c(c(0,0,1), vector.product(direction[i,], c(0,0,1))))))
  src <- sweep(-DSA * direction, 2, PTV.G,FUN = "+")
  
  src.extrem.pt.dist <- max(apply(src,1, function(sc)  
    max(.fnorm(sweep(t((vol$xyz.from.ijk %*% vol$cube.idx)[1:3,]), 2, sc)))))
  dalpha <- 180 * min(vol$dxyz) / src.extrem.pt.dist / pi
  alpha <- alpha + dalpha
  
  density <- vol.copy(vol)
  density$vol3D.data[density$vol3D.dat <= -900]  <- -1000
  density$vol3D.data[]  <- (density$vol3D.data[] / 1000 + 1)
  density$max.pixel <- (density$max.pixel / 1000 + 1)
  density$min.pixel <- (density$min.pixel / 1000 + 1)
  
  D.list <- lapply( 1:beam.nb, function(theta.idx){
    
    fan <- fan.beam (alpha =alpha,dalpha, origin = as.numeric(src[theta.idx,]),
                     orientation =  as.numeric(orientation[theta.idx,]),
                     ref.pseudo = bin.PTV.d$ref.pseudo)
    vox <- fan.to.voxel(bin.PTV.d,fan,restrict = T )
    ray.index <- sort( unique(vox$ray.index))
    fan$xyz <- fan$xyz[ray.index, ]
    fan$local.coord <- fan$local.coord[ray.index, ]
    r <- as.numeric(DSA*sqrt((fan$xyz %*% fan$orientation[1:3])^2 + (fan$xyz %*% fan$orientation[4:6])^2))
    flou <- 1 - pnorm(r,tan(alpha*pi/180) * DSA * 0.9 , 1)
    voxel <- .fan.to.voxel.with.att(density,fan, att = TRUE)
    voxel$att  <-  pnorm(voxel$att, 2, 5) * exp(-voxel$att / 150)
    vol.idx <- unique(voxel$vol.index)
    D_ <- vol.copy(vol, modality = "rtdose")
    D_$vol3D.data[] <- 0
    D_$vol3D.data <- .mean_voxC(as.numeric(D_$vol3D.data), vol.idx - 1, 
                                voxel$vol.index - 1,voxel$att, flou[voxel$ray.index])
    D_$vol3D.data  <- array(D_$vol3D.data ,dim = vol$n.ijk)
    D_$vol3D.data  <- D_$vol3D.data / max(D_$vol3D.data,na.rm = TRUE)
    D_$max.pixel <- 1
    D_$min.pixel <- 0
    D_
  })
  
  D <- vol.copy(vol,modality = "rtdose" )
  D$vol3D.data[] <- 0
  D$min.pixel <- 0
  D$max.pixel <- 0
  for (i in 1:beam.nb) D <- vol.sum(D,D.list[[i]], alias = alias, description = description)
  
  D$vol3D.data <- D.max*D$vol3D.data / nesting.cube(D, pt.min = pt.min,
                                                    pt.max =  pt.max)$max.pixel
  D$max.pixel  <- max(D$vol3D.data,na.rm = TRUE)
  D$min.pixel  <- min(D$vol3D.data,na.rm = TRUE)
  return(D)
}