R/transform3DIMS.R

In NHPatterson/Reg3DIMS: Automated registration of sequential 2D IMS for 3D reconstruction in Cardinal

#' @export
#' @title Transfom 2D IMS data into 3D reconstruction
#' @description Generates a name ordered list of MSImagingExperiment data that has been transformed into 3D
#' @param ims Cardinal data of class: MSImagingExperiment
#' @param tforms named list of NiftyReg transformations per 2D section, generated by reg3DIMS function
#' @return alphanumerically ordered named list of MSImagingExperiment's transformed to their 3D coordinates

transform3DIMS <- function(ims, tforms){
  ims_data <- list()
  nsignals <- nrow(ims)[[1]]
  for(im in names(tforms)){
    cat('Transforming :', im,'\n')
    tslice <- slice(ims[,run(ims) == im],feature=1)
    
    if(is.null(dim(tforms[[im]]$nifty$image))){
      cat('mid point image','\n')
      mpdata <- ims[,run(ims) == im]
      xpad_left <- floor((canvas_xy[1] - max(coord(mpdata)$x)) / 2)
      ypad_top  <- floor((canvas_xy[2] - max(coord(mpdata)$y)) / 2)
      coord(mpdata)$x <- coord(mpdata)$x + xpad_left
      coord(mpdata)$y <- coord(mpdata)$y + ypad_top
      
      ims_data[[im]] <- mpdata
    }else{
      canvas_xy <- dim(tforms[[im]]$nifty$image)
      
      #calculate padding
      xpad_left <- floor((canvas_xy[1] - nrow(tslice)) / 2)
      ypad_top <- floor((canvas_xy[2] - ncol(tslice)) / 2)
      
      #binary slice
      tslice[!is.na(tslice)] <- 1
      tslice[is.na(tslice)] <- 0
      binaryslice <- array(rep(0, canvas_xy[1]*canvas_xy[2]), dim=c(canvas_xy[1],canvas_xy[2]))
      
      binaryslice[xpad_left:(xpad_left+nrow(tslice) -1), ypad_top:(ypad_top+ncol(tslice) -1)] <- tslice

      #empty array
      mzslice <- array(rep(0, canvas_xy[1]*canvas_xy[2]*nrow(ims)[[1]]), dim=c(canvas_xy[1],canvas_xy[2],nsignals))
      cat('Extracing IMS data cube','\n')
      
      #fill array by cardinal slice
      mzs <- Cardinal::mz(ims)
      mzslice[xpad_left:(xpad_left+nrow(tslice) -1), ypad_top:(ypad_top+ncol(tslice) -1),] <- slice(ims[,run(ims) == im], mz = mzs )
      
      mzslice[is.na(mzslice)] <- 0
      
      
      cat('Transforming IMS data cube','\n')
      for(signal in seq_len(dim(mzslice)[3])){
        if(!is.null(tforms[[im]]$init)){
          mzslice[,,signal] <- applyTransform(tforms[[im]]$init, mzslice[,,signal],interpolation = 0)
        }
        mzslice[,,signal] <- applyTransform(forward(tforms[[im]]$nifty), mzslice[,,signal],interpolation = 0)
        
      }
      
      if(!is.null(tforms[[im]]$init)){
        binaryslice_tformed <- applyTransform(tforms[[im]]$init, binaryslice,interpolation = 0)
        binaryslice_tformed <- applyTransform(forward(tforms[[im]]$nifty), binaryslice_tformed,interpolation = 0)
        
      }else{
        binaryslice_tformed <- applyTransform(forward(tforms[[im]]$nifty), binaryslice,interpolation = 0)
        
      }
      
      
      xy <- data.frame(x=as.vector(row(binaryslice_tformed)), y=as.vector(col(binaryslice_tformed)))
      retained_idx <- which(as.vector(binaryslice_tformed) == 1)
      
      xy <- xy[retained_idx,]
      
      cat('Rebuliding transformed IMS data','\n')
      
      for(signal in seq_len(dim(mzslice)[3])){
        intdatabuffer <- as.matrix(as.vector(mzslice[,,signal])[retained_idx])
        if(signal > 1){
          intdata <- cbind(intdata, intdatabuffer)
        }else{
          intdata <- intdatabuffer
        }

      }
      

      ims_data[[im]] <- MSImagingExperiment(imageData = t(as.matrix(intdata)),
                          featureData = MassDataFrame(mz=Cardinal::mz(ims)),
                          pixelData = PositionDataFrame(coord = xy, run = im))
      
      }
    
  }
  ims_data <- ims_data[order(names(ims_data))]
  return(ims_data)
}