Wishlist: add mgh to nii conversion

[dipterix]

Hi @dfsp-spirit , I wonder how hard it would be to write a mgh/mgz to nii converter? I tried to implement by myself, but I don't know where to get rest of the header information or if the conversion is generally applicable. Can I steal some knowledge from you?

Thanks!

m44_to_quaternion <- function(m) {
  m00 <- m[1, 1]
  m01 <- m[1, 2]
  m02 <- m[1, 3]
  m10 <- m[2, 1]
  m11 <- m[2, 2]
  m12 <- m[2, 3]
  m20 <- m[3, 1]
  m21 <- m[3, 2]
  m22 <- m[3, 3]

  tr <- m00 + m11 + m22 + 1.0

  # https://github.com/NIFTI-Imaging/nifti_clib/blob/d05a71ecdc0b811a509b77c3e82fd0674e901807/niftilib/nifti1_io.c
  if( tr > 0.5 ) {
    S <- sqrt( tr ) * 2.0
    qw <- 0.25 * S
    qx <- (m21 - m12) / S
    qy <- (m02 - m20) / S
    qz <- (m10 - m01) / S
  } else {
    Sx <- sqrt(1.0 + m00 - m11 - m22) * 2.0 # S = 4 * qx
    Sy <- sqrt(1.0 + m11 - m00 - m22) * 2.0 # S = 4 * qy
    Sz <- sqrt(1.0 + m22 - m00 - m11) * 2.0 # S = 4 * qz

    if( Sx > 2.0 ) {
      qw <- (m21 - m12) / Sx
      qx <- 0.25 * Sx
      qy <- (m01 + m10) / Sx
      qz <- (m02 + m20) / Sx
    } else if ( Sy > 2.0 ) {
      qw <- (m02 - m20) / Sy
      qx <- (m01 + m10) / Sy
      qy <- 0.25 * Sy
      qz <- (m12 + m21) / Sy
    } else {
      qw <- (m10 - m01) / Sz
      qx <- (m02 + m20) / Sz
      qy <- (m12 + m21) / Sz
      qz <- 0.25 * Sz
    }

    if( tr < 0.0 ) {
      qx <- -qx
      qy <- -qy
      qz <- -qz
    }

  }
  return(c(qw, qx, qy, qz))
}

mgh_to_nii <- function(from, to) {

  # from <- "/Users/dipterix/Dropbox (PennNeurosurgery)/RAVE/Samples/raw/PAV006/rave-imaging/fs/mri/aparc+aseg.mgz"
  # to <- tempfile(fileext = ".nii")

  mgh <- freesurferformats::read.fs.mgh(from, with_header = TRUE, drop_empty_dims = FALSE)

  # generate nii1 header
  header <- freesurferformats:::ni1header.for.data(mgh$data)
  header$dim <- c(3L, dim(mgh$data)[1:3], 1L, 1L, 1L, 1L)
  header$intent_p1 <- 0
  header$intent_p2 <- 0
  header$intent_p3 <- 0
  header$intent_code <- 0L
  # header$datatype <- 8L
  # header$bitpix <- 32L
  header$slice_start <- 0L
  header$pix_dim <- c(-1, 1, 1, 1, 0, 1, 1, 1)
  # header$vox_offset <- 352L
  # header$scl_slope
  # header$scl_inter
  # header$slice_end
  # header$slice_code
  header$xyzt_units <- 10L
  # header$cal_max
  # header$cal_min
  # header$slice_duration <- 0L
  # header$toffset
  # header$glmax
  header$qform_code <- 1L
  mat <- mgh$header$vox2ras_matrix
  # a = 0.5  * sqrt(1+R11+R22+R33)    (not stored)
  # b = 0.25 * (R32-R23) / a       => quatern_b
  # c = 0.25 * (R13-R31) / a       => quatern_c
  # d = 0.25 * (R21-R12) / a       => quatern_d

  # extract rotation
  rot <- mat
  rot[, 4] <- c(0,0,0,1)
  rot <- t(t(rot) / sqrt(colSums(rot^2)))
  if( det(rot) < 0 ) {
    # r13 = -r13 ; r23 = -r23 ; r33 = -r33 ;
    rot[,3] <- -rot[,3]
  }
  quatern <- m44_to_quaternion(rot)
  header$quatern_b <- quatern[2]
  header$quatern_c <- quatern[3]
  header$quatern_d <- quatern[4]

  header$sform_code <- 1L
  header$qoffset_x <- mat[1, 4]
  header$qoffset_y <- mat[2, 4]
  header$qoffset_z <- mat[3, 4]
  header$srow_x <- mat[1, ]
  header$srow_y <- mat[2, ]
  header$srow_z <- mat[3, ]

  freesurferformats::write.nifti1(
    filepath = to,
    niidata = mgh$data,
    niiheader = header
  )

}

[dfsp-spirit]

Converting MGH/MGZ brain volumes (3D/4D) to NIFTI brain volumes is definitely possible, and you can get most of the required information from a proper MGH file. There are some limitations though.

Note: Of course, some software libraries allow people to store data without proper header (basically the MGH is just a 4D array then), or some people store 1D-data (per-vertex overlays) in MGH files instead of using curv format, which also have no useful header then, of course. But if the MGH is a proper file representing a volume image, the information is mostly there. You can tell whether that's the case from the is_ras_valid field in the MGH file.

MGH has, for example, no INTENT field (see NIFTI1 standard), so your converter would have to assume that the source MGH/MGZ is a brain volume (people could store whatever in the file, but you cannot know what they meant) and set some NIFTI header fields to zeros that are not useful for brain volumes. If the MGH file contains data with a different meaning (like p-values for each voxel, so a stats map), you cannot really know that, and so you cannot properly set the INTENT field in general. You have to make some assumptions. (Or you need a way more complex function that takes as an additional argument an INTENT, and fills the NIFTI header information based on the data and that intent information.)

You can see the inverse operation to what you want to do here in nifti_to_mgh.R.

Also note that various NIFTI header fields are simply unused (see linked NIFTI 1 standard document, Chapter 7).

I thought I had already written a write_fs_vol() function that could easily save as nifti, but apparently not. I think what freesurferformats needs is a function freesurferformats:::ni1header.for.mgh(mgh), in addition to the existing freesurferformats:::ni1header.for.data(mgh$data). Then you would basically have your converter.

It's been a while since I messed with the NIFTI standard, and I would have to read up on the details of the header fields myself I guess.

[dipterix]

freesurferformats:::ni1header.for.mgh sounds like a good a idea.

[dfsp-spirit]

So it is fine with you if we assume that the file contains a 3D/4D brain volume, right? Then I could have a look.

[dipterix]

Right. Most likely when people read a mgh file, they are reading mri/ folder, which contains a lot of T1, aparc aseg, ....
Even assuming the file containing a 3D/4D brain volume could be very beneficial.

I guess as long as you mention in the help document that this is only intended for 3D/4D brain volumes, people should respect it.

Thanks a lot : )

[dfsp-spirit]

Yes, what people want in the end, in the case of a 3D/4D volume, is the correct RAS information to be extracted from the NII and stored in the respective header fields of the MGH.