R/nifti_io.R

In neuroim2: Data Structures for Brain Imaging Data

#' @include common.R
{}
#' @include binary_io.R
{}
#' @include nifti_extensions.R
{}
#' @keywords internal
#' @noRd
.checkDimensions <- function(dimvec) {
  if (!all(dimvec >= 0)) {
    cli::cli_abort("Illegal dimension vector in header: {paste(dimvec, collapse=' x ')}.")
  }
}


#' @keywords internal
#' @noRd
write_nifti_vector <- function(vec, file_name, data_type="FLOAT") {
  if (length(dim(vec)) != 4) {
    cli::cli_abort("Input vector must be 4-dimensional, not {length(dim(vec))}D.")
  }
	ext <- .make_nifti_volume_labels_extension(volume_labels(vec))
	hdr <- as_nifti_header(vec, file_name = file_name, data_type = data_type,
                         extensions = ext)

	conn <- if (substr(file_name, nchar(file_name)-2, nchar(file_name)) == ".gz") {
				gzfile(file_name, open="wb")
			} else {
				file(file_name, open="wb")
			}

	write_nifti_header(hdr, conn, close=FALSE)
	writer <- BinaryWriter(conn, hdr$vox_offset, data_type, hdr$bitpix/8, .Platform$endian)

	NVOLS <- dim(vec)[4]

	for (i in 1:NVOLS) {
		write_elements(writer, as.numeric(vec[[i]]))
	}
	close(writer)
}

#' @keywords internal
#' @noRd
write_nifti_hyper_vector <- function(vec, file_name, data_type="FLOAT") {
  if (!inherits(vec, "NeuroHyperVec")) {
    cli::cli_abort("{.arg vec} must be a {.cls NeuroHyperVec} object.")
  }
  if (length(dim(vec)) != 5) {
    cli::cli_abort("Input hyper-vector must be 5-dimensional, not {length(dim(vec))}D.")
  }

	hdr <- createNIfTIHeader(oneFile=TRUE, file_name=file_name)
	hdr$file_name <- file_name
	hdr$endian <- .Platform$endian
	hdr$datatype <- .getDataCode(data_type)
	hdr$data_storage <- .getDataStorage(hdr$datatype)
	hdr$bitpix <- .getDataSize(data_type) * 8
	hdr$dimensions <- c(length(dim(vec)), dim(vec))
	N <- 8 - length(hdr$dimensions)
	hdr$dimensions <- c(hdr$dimensions, rep(1, N))
	hdr$num_dimensions <- length(dim(vec))

	hdr$pixdim <- c(0, spacing(vec), rep(1, 4))
	hdr$scl_intercept <- 0
	hdr$scl_slope <- 1

	tmat <- trans(vec)
	hdr$qoffset <- tmat[1:3, 4]
	hdr$qform <- tmat
	hdr$sform <- tmat
	quat1 <- matrixToQuatern(tmat)
	hdr$quaternion <- quat1$quaternion
	hdr$qfac <- quat1$qfac
	hdr$pixdim[1] <- hdr$qfac

	hdr$extensions <- new("NiftiExtensionList")
	hdr$vox_offset <- 348 + total_extension_size(hdr$extensions)

	conn <- if (substr(file_name, nchar(file_name)-2, nchar(file_name)) == ".gz") {
		gzfile(file_name, open="wb")
	} else {
		file(file_name, open="wb")
	}

	write_nifti_header(hdr, conn, close=FALSE)
	writer <- BinaryWriter(conn, hdr$vox_offset, data_type, hdr$bitpix/8, .Platform$endian)
	write_elements(writer, as.numeric(dense_array_5d(vec)))
	close(writer)
}

#' @keywords internal
#' @noRd
write_nifti_volume <- function(vol, file_name, data_type="FLOAT") {
  if (length(dim(vol)) != 3) {
    cli::cli_abort("Input volume must be 3-dimensional, not {length(dim(vol))}D.")
  }
	hdr <- as_nifti_header(vol, file_name=file_name, data_type=data_type)

	is_gzipped <- endsWith(file_name, ".gz")
	conn <- if (is_gzipped) {
		gzfile(file_name, open="wb")
	} else {
		file(file_name, open="wb")
	}

	write_nifti_header(hdr, conn, close=FALSE)
	writer <- BinaryWriter(conn, hdr$vox_offset, data_type, hdr$bitpix/8, .Platform$endian)
	write_elements(writer, as.numeric(vol))
	close(writer)
}


#' Construct a Minimal NIfTI-1 Header from a NeuroVol
#'
#' @description
#' Given a \code{\link[neuroim2]{NeuroVol}} object (or similar), this function
#' builds a basic NIfTI-1 header structure, populating essential fields such as
#' \code{dim}, \code{pixdim}, \code{datatype}, the affine transform, and the
#' quaternion parameters.
#'
#' @details
#' This is a convenience function that calls \code{\link{createNIfTIHeader}}
#' first, then updates the fields (dimensions, \code{pixdim}, orientation, etc.)
#' based on the \code{vol} argument. The voxel offset is set to 352 bytes (or
#' larger if extensions are provided), and the quaternion is derived from the
#' transform matrix via \code{\link{matrixToQuatern}}.
#'
#' Note: This function primarily sets up a minimal header suitable for writing
#' standard single-file NIfTI-1. If you need a more comprehensive or advanced
#' usage, consider manually editing the returned list.
#'
#' @param vol A \code{\link[neuroim2]{NeuroVol}} (or 3D array-like) specifying
#'   dimensions, spacing, and affine transform.
#' @param file_name A character string for the file name (used within the header
#'   but not necessarily to write data).
#' @param oneFile Logical; if \code{TRUE}, sets the NIfTI magic to \code{"n+1"},
#'   implying a single-file format (\code{.nii}). If \code{FALSE}, uses
#'   \code{"ni1"} (header+image).
#' @param data_type Character specifying the data representation, e.g. \code{"FLOAT"},
#'   \code{"DOUBLE"}. The internal code picks an integer NIfTI code.
#' @param extensions Optional \code{\link{NiftiExtensionList-class}} object or list
#'   of \code{\link{NiftiExtension-class}} objects to include in the header.
#'
#' @return A \code{list} representing the NIfTI-1 header fields, containing
#'   elements like \code{dimensions}, \code{pixdim}, \code{datatype}, \code{qform},
#'   \code{quaternion}, \code{qfac}, \code{extensions}, etc. This can be passed to
#'   other functions that write or manipulate the header.
#'
#' @seealso
#' \code{\link{createNIfTIHeader}} for the base constructor of an empty NIfTI header.
#' \code{\link{NiftiExtension}} for creating extensions.
#'
#' @export
as_nifti_header <- function(vol, file_name, oneFile=TRUE, data_type="FLOAT",
                           extensions = NULL) {
		hd <- createNIfTIHeader(oneFile=oneFile, file_name=file_name)
		hd$file_name <- file_name
		hd$endian <- .Platform$endian
		hd$datatype <- .getDataCode(data_type)
		hd$data_storage <- .getDataStorage(hd$datatype)
		hd$bitpix <- .getDataSize(data_type) * 8
		hd$dimensions <- c(length(dim(vol)), dim(vol))
		N <- 8 - length(hd$dimensions)
		hd$dimensions <- c(hd$dimensions,  rep(1, N))
		hd$num_dimensions <- length(dim(vol))

		### only encodes pixdim for three dimensions
		hd$pixdim <- c(0, spacing(vol), rep(0,4))

		hd$scl_intercept <- 0
		hd$scl_slope <- 1

		tmat <- trans(vol)

		# Derive qoffset from the transform matrix translation column
		# to guarantee consistency between qoffset and sform
		hd$qoffset <- tmat[1:3, 4]

		hd$qform <- tmat
		hd$sform <- tmat

		quat1 <- matrixToQuatern(tmat)
		hd$quaternion <- quat1$quaternion
		hd$qfac <- quat1$qfac
		hd$pixdim[1] <- hd$qfac

		# Handle extensions
		if (!is.null(extensions)) {
		  if (is.list(extensions) && !is(extensions, "NiftiExtensionList")) {
		    # Convert list of NiftiExtension objects to NiftiExtensionList
		    extensions <- new("NiftiExtensionList", extensions)
		  }
		  hd$extensions <- extensions
		} else {
		  hd$extensions <- new("NiftiExtensionList")
		}

		# Calculate vox_offset: 348 (header) + extension bytes
		# total_extension_size includes the 4-byte extender field
		hd$vox_offset <- 348 + total_extension_size(hd$extensions)

		hd

}


#' Create an Empty NIfTI-1 Header List
#'
#' @description
#' Initializes a list of fields following the NIfTI-1 specification with default
#' or placeholder values. Users typically call this internally via
#' \code{\link{as_nifti_header}} rather than using directly.
#'
#' @details
#' This function sets up the skeleton of a NIfTI-1 header, including fields for
#' \code{diminfo}, \code{pixdim}, \code{qform_code}, \code{magic}, etc. Most fields
#' are initialized to zero, empty characters, or standard placeholders. The
#' \code{oneFile} argument controls whether \code{"n+1"} or \code{"ni1"} is used
#' for the \code{magic} field.
#'
#' @param oneFile Logical; if \code{TRUE}, \code{magic} is set to \code{"n+1"}
#'   indicating a single-file (.nii) approach. Otherwise set to \code{"ni1"}.
#' @param file_name Optional character string to store in the header, usually
#'   referencing the intended output file name.
#'
#' @return A named \code{list} containing approximately 30 fields that comprise
#'   the NIfTI-1 header structure. Many of these are placeholders until filled by
#'   downstream usage.
#'
#' @seealso
#' \code{\link{as_nifti_header}} for populating the returned header with actual
#' data from a \code{NeuroVol}.
#'
#' @export
createNIfTIHeader <- function(oneFile=TRUE, file_name=NULL) {
	header <- list()
	header$file_type <- "NIfTI"
	header$encoding <- "binary"
	header$version <- "1"


	header$file_name <- file_name
	header$endian <- .Platform$endian

	header$diminfo <- 0
	header$dimensions <- NULL
	header$num_dimensions <- NULL

	header$intent1 <-  0
	header$intent2 <-  0
	header$intent3 <-  0


	header$intent_code <-  0
	header$datatype <- NULL
	header$data_storage <- NULL
	header$bitpix <- NULL
	header$slice_start <- 0
	header$pixdim <-  NULL

	header$qfac <- -1


	header$vox_offset <- 0
	header$scl_slope <- 0
	header$scl_intercept <- 0
	header$slice_end <- 0
	header$slice_code <- 0


	header$xyzt_units <- 2
	header$cal_max <- 0
	header$cal_min <- 0

	header$slice_duration <- 0
	header$toffset <- 0

	header$glmax <- 0
	header$glmin <- 0

	header$description <- character(80)
	header$auxfile <- character(24)

	header$qform_code <- 1
	header$sform_code <- 1
	header$quaternion <- NULL

	header$qoffset <- NULL
	header$qform <- NULL

	header$sform <- NULL
	header$intent_name <- character(16)
	header$magic <- "n+1"

	header$onefile <- oneFile
	if (oneFile) {
		header$magic <- "n+1"
	} else {
		header$magic <- "ni1"
	}

	header$version <- 1
	header

}


#' Read NIfTI Header
#'
#' @description
#' Reads a NIfTI-1 header from file, including any extensions present.
#'
#' @param fname Character string specifying the path to the NIfTI file
#'   (.nii, .nii.gz, or .hdr).
#' @param read_extensions Logical; if TRUE (default), reads any extensions
#'   present in the file. Set to FALSE to skip extension reading.
#'
#' @return A list containing header fields including:
#'   \itemize{
#'     \item Standard NIfTI header fields (dimensions, pixdim, transforms, etc.)
#'     \item \code{extensions}: A \code{\link{NiftiExtensionList-class}} object
#'       containing any extensions found (empty if none present or read_extensions=FALSE)
#'   }
#'
#' @keywords internal
#' @noRd
read_nifti_header <- function(fname, read_extensions = TRUE) {
	header <- list()
	header$file_type <- "NIfTI"
	header$encoding <- "binary"
	header$version <- "1"

	conn <- NULL

	if (.isExtension(fname, ".nii") || .isExtension(fname, ".hdr")) {
		conn <- file(fname, open="rb")
	} else if (.isExtension(fname, ".nii.gz")) {
		conn <- gzfile(fname, open="rb")
		header$encoding <- "gzip"
	} else {
		stop(paste("illegal NIFTI header name", fname))
	}

	# Ensure connection is closed on exit
	on.exit(close(conn), add = TRUE)

	endian <- .getEndian(conn)

	header$file_name <- fname
	header$endian <- endian

	readBin(conn, what=integer(), n=10+18+4+2+1, size=1)

	header$diminfo <- readBin(conn, what=integer(), n=1, size=1)
	header$dimensions <- readBin(conn, integer(), n=8, size=2, endian=endian)

	header$dimensions[header$dimensions == 0] <- 1

	.checkDimensions(header$dimensions)

	header$num_dimensions <- header$dimensions[1]

	if (header$num_dimensions == 5 && header$dimensions[5] == 1) {
	  warning("read_nifti_header: collapsing fourth dimension because it has length = 1")
	  newdim <- header$dimensions[c(1:4, 6)]
	  header$num_dimensions <- 4
	  header$dimensions <- c(newdim, rep(1,3))
	  header$dimensions[1] <- 4
	}

	header$intent1 <-  readBin(conn, double(), n=1, size=4, endian=endian)
	header$intent2 <-  readBin(conn, double(), n=1, size=4, endian=endian)
	header$intent3 <-  readBin(conn, double(), n=1, size=4, endian=endian)


	header$intent_code <-  readBin(conn, integer(), n=1, size=2, endian=endian)
	header$datatype <- readBin(conn, integer(), n=1, size=2, endian=endian)
	header$data_storage <- .getDataStorage(header$datatype)
	header$bitpix <- readBin(conn, integer(), n=1, size=2, endian=endian)
	header$slice_start <- readBin(conn, integer(), n=1, size=2, endian=endian)
	header$pixdim <-  readBin(conn, double(), n=8, size=4, endian=endian)

	header$qfac = header$pixdim[1]

	if (header$qfac == 0) {
		header$qfac = 1
	}

	header$vox_offset <- readBin(conn, double(), n=1, size=4, endian=endian)
	header$scl_slope <- readBin(conn, double(), n=1, size=4, endian=endian)
	header$scl_intercept <- readBin(conn, double(), n=1, size=4, endian=endian)
	header$slice_end <- readBin(conn, integer(), n=1, size=2, endian=endian)
	header$slice_code <-  readBin(conn, integer(), n=1, size=1, endian=endian)


	header$xyzt_units <- readBin(conn, integer(), n=1, size=1, endian=endian)
	header$cal_max <- readBin(conn, double(), n=1, size=4, endian=endian)
	header$cal_min <- readBin(conn, double(), n=1, size=4, endian=endian)

	header$slice_duration <- readBin(conn, double(), n=1, size=4, endian=endian)
	header$toffset <- readBin(conn, double(), n=1, size=4, endian=endian)

	header$glmax <- readBin(conn, integer(), n=1, size=4, endian=endian) # unused glmax, glmin
	header$glmin <- readBin(conn, integer(), n=1, size=4, endian=endian) # unused glmax, glmin

	header$description <- readBin(conn, integer(), n=80, size=1, endian=endian)
	header$auxfile <- readBin(conn, integer(), n=24, size=1, endian=endian)

	header$qform_code <- readBin(conn, integer(), n=1, size=2, endian=endian)
	header$sform_code <- readBin(conn, integer(), n=1, size=2, endian=endian)
	header$quaternion <- readBin(conn, double(), n=3, size=4, endian=endian)

	header$qoffset <- readBin(conn, double(), n=3, size=4, endian=endian)
	header$qform <- quaternToMatrix(header$quaternion, header$qoffset, header$pixdim[2:4], header$qfac)

	sform  <- readBin(conn, double(), n=12, size=4, endian=endian)
	header$sform <- rbind(matrix(sform,3,4, byrow=T), c(0,0,0,1))
	header$intent_name <- readBin(conn, character(), n=16, size=1, endian=endian)
	header$magic <- suppressWarnings(readChar(conn, nchars=4))

	header$onefile <- F
	if (substr(header$magic,2,2) == "+") {
		header$onefile <- T
	}

	header$version <- substr(header$magic,3,3)

	# Read extensions if requested and if this is a single-file NIfTI
	# Extensions are only valid in single-file format (.nii) and when vox_offset > 352
	if (read_extensions && header$onefile && header$vox_offset > 352) {
		# Connection is at byte 348 after reading magic (344) + 4 bytes
		header$extensions <- read_nifti_extensions(conn, header$vox_offset, endian)
	} else {
		header$extensions <- new("NiftiExtensionList")
	}

	header

}


#' Write NIfTI Header
#'
#' @description
#' Writes a NIfTI-1 header to a connection, including any extensions.
#'
#' @param niftiInfo A list containing NIfTI header fields.
#' @param conn An open binary connection to write to.
#' @param close Logical; if TRUE (default), close the connection after writing.
#'
#' @return The connection (invisibly).
#'
#' @details
#' The function writes the 348-byte standard header, then extensions (if present),
#' then pads to vox_offset. The vox_offset in niftiInfo should already account
#' for extension size.
#'
#' @keywords internal
#' @noRd
write_nifti_header <- function(niftiInfo, conn, close=TRUE) {
	endian <- niftiInfo$endian


	writeBin(as.integer(348), conn, 4, endian)
	writeBin(integer(34),conn,1,endian)
	writeChar("r", conn,1,eos=NULL)
	writeBin(as.integer(niftiInfo$diminfo), conn, size=1, endian) #diminfo, not supported currently -- write zero
	#writeBin(as.integer(niftiInfo$num_dimensions), conn, 2, endian)         #num dimensions

	stopifnot(length(niftiInfo$dimensions) == 8)
	stopifnot(length(niftiInfo$pixdim) == 8)

	writeBin(as.integer(niftiInfo$dimensions), conn, 2, endian)   #dimension vector
	writeBin(as.double(niftiInfo$intent1), conn, 4, endian)       #intent1
	writeBin(as.double(niftiInfo$intent2), conn, 4, endian)       #intent2
	writeBin(as.double(niftiInfo$intent3), conn, 4, endian)       #intent3
	writeBin(as.integer(niftiInfo$intent_code), conn, 2, endian)   #intent code
	writeBin(as.integer(niftiInfo$datatype),conn, 2, endian)      #datatype
	writeBin(as.integer(niftiInfo$bitpix),conn, 2, endian)        #bits per pixel
	writeBin(as.integer(niftiInfo$slice_start),conn, 2, endian)    #slice start
	writeBin(as.double(niftiInfo$pixdim), conn, 4, endian)        #pix dim
	writeBin(as.double(niftiInfo$vox_offset), conn, 4, endian)     #voxel offset
	writeBin(as.double(niftiInfo$scl_slope), conn, 4, endian)      #slope
	writeBin(as.double(niftiInfo$scl_intercept), conn, 4, endian)  #intercept
	writeBin(as.integer(niftiInfo$slice_end), conn, 2, endian)     #slice end
	writeBin(as.integer(niftiInfo$slice_code), conn, 1, endian)    #slice code
	writeBin(as.integer(niftiInfo$xyzt_units), conn, 1, endian)    #xyzt units
	writeBin(as.double(niftiInfo$cal_max), conn, 4, endian)        #cal max
	writeBin(as.double(niftiInfo$cal_min), conn, 4, endian)        #cal min
	writeBin(as.double(niftiInfo$slice_duration), conn, 4, endian) #slice duration
	writeBin(as.double(niftiInfo$toffset), conn, 4, endian)       #t offset
	writeBin(as.integer(niftiInfo$glmax), conn, 4, endian)        #glmax
	writeBin(as.integer(niftiInfo$glmin), conn, 4, endian)        #glmin
	writeBin(as.integer(niftiInfo$description), conn, 1, endian)  #description
	writeBin(as.integer(niftiInfo$auxfile), conn, 1, endian)      #aux_file
	writeBin(as.integer(niftiInfo$qform_code), conn, 2, endian)    #qform code
	writeBin(as.integer(niftiInfo$sform_code), conn, 2, endian)    #sform code

	writeBin(as.double(niftiInfo$quaternion), conn, 4, endian)    #quaternion
	writeBin(as.double(niftiInfo$qoffset), conn, 4, endian)       #qoffset
	writeBin(as.double(t(niftiInfo$sform[1:3,])), conn, 4, endian) #sform
	writeBin(as.integer(niftiInfo$intent_name), conn, 1, endian)    #intent_name
	writeChar(niftiInfo$magic, conn)                               #magic

	# We're now at byte 348
	# Write extensions if present
	extensions <- niftiInfo$extensions
	write_nifti_extensions(conn, extensions, endian)

	# Pad remaining bytes to vox_offset if needed
	loc <- seek(conn)
	offset <- niftiInfo$vox_offset

	nbytes <- offset - loc
	if (nbytes > 0) {
		writeBin(raw(nbytes), conn)
	}

	if (close) {
		close(conn)
	}

	invisible(conn)
}