Cannot resample very high-res surfaces (due to encoding of written GIFTI files)

[smeisler]

Hello,

Crosspost of mandymejia/BayesfMRI#16

All relevant surface fmriprep outputs for a single subject, as well as the 32k surface fsLR template from templateflow, can be found here.

The goal is to resample the high-density subject surface (sub-047EPKL011005_ses-1_hemi-L_midthickness.surf.gii) to the same resolution as the template (tpl-fsLR_den-32k_hemi-L_midthickness.surf.gii), which matches the cifti (e.g. sub-047EPKL011005_ses-2_task-read_run-1_space-fsLR_den-91k_bold.dtseries.nii).

I also included things like confounds and events, in case these data need to be tested further on BayesfMRI. These data are public (openneuro dataset ds003126), so I have no concerns with sharing derivatives.

Command

test <- resample_surf(surf=fname_gifti_left, resamp_res = 32492, hemisphere='left')

Error

Error in read_surf(surf = surf, expected_hemisphere = expected_hemisphere, : The object could not be converted into a surface.
Traceback:

1. resample_surf(surf = fname_gifti_left, resamp_res = 32492, hemisphere = "left")
2. make_surf(gii_post, hemisphere)
3. read_surf(surf = surf, expected_hemisphere = expected_hemisphere, 
 .     resamp_res = NULL)
4. stop("The object could not be converted into a surface.")

And (I cut out most of the error message because it is very big and repetitive, the last part is below):
triangle 317355, vertex 0Vertex used twice in one triangle: triangle 317356, vertex 0Vertex used twice in one triangle: triangle 317356, vertex 0Vertex used twice in one triangle: triangle 317357, vertex 0Vertex used twice in one triangle: triangle 317357, vertex 0...

Please let me know if I can provide more info or help in any way.

Thanks,
Steven

[damondpham]

Thank you Steven, this is very helpful!

You mentioned that you have sphere files in the same resolution as the surface, sub-047EPKL011005_ses-1_hemi-L_midthickness.surf.gii. Could you also attach those here?

Part of the problem seems to be that the Connectome Workbench command "-surface-resample", which ciftiTools depends on, does not generate surfaces with the same resolution as this surface. For example, resampling a default surface at a requested resolution of 158681 yields a surface with a resolution of 158762.

[smeisler]

Sure! Those files are at https://drive.google.com/drive/folders/13oLnnjtnD0fpMN1InQZEygSShNFLjXyM?usp=sharing

[damondpham]

Thanks! It looks like the spheres are actually in 32k, whereas the surface is in ~150k?

> surfL_fname <- "sub-047EPKL011005_ses-1_hemi-L_midthickness.surf.gii"
> sphereL_fname <- "tpl-fsLR_hemi-L_den-32k_sphere.surf.gii"
> read_surf(surfL_fname)
Vertices:    158681 
Faces:       317358 
Hemisphere:  left 
Warning message:
In rbind(c(1.045321, NA, NA, 0.015827, NA, NA, 0.068957, NA, NA,  :
  number of columns of result is not a multiple of vector length (arg 2)
> read_surf(sphereL_fname)
Vertices:    32492 
Faces:       64980 
Hemisphere:  left 

(The warning is caused by metadata in the surface GIFTI, and occurs with read_gifti(surfL_fname), so it's a separate issue I'll look into)

[smeisler]

Sorry for the misunderstanding, I meant originally that the spheres match the template surface and func cifti, not the high-res subject surface.

[damondpham]

I see. Either way, it looks like the Connectome Workbench command "-surface-resample" cannot resample an input surface that is very-high resolution. This example causes the same problem:

surfDefL <- load_surf(resamp_res=150000) # create a 150k input surface from the default geometry
resample_surf(surfDefL, 1000) # doesn't work

Note that a surface can be resampled to a very-high resolution; the error seems to occur when a surface is resampled from a very-high resolution to any other resolution, larger or smaller.

I will have to look into this more, but it seems like for now, there isn't a way that ciftiTools can work with surfaces that are very-high resolution, unless a lower resolution version already exists.

[smeisler]

I wonder if it has to do with the temporary file location, I found this at the top of the error message after scrolling past the zero triangle warnings:

While running:
/rdma/vast-rdma/vast/gablab/smeisler/workbench/bin_rh_linux64/../exe_rh_linux64/wb_command -surface-resample /tmp/Rtmp49kVVm/to_resample.surf.gii /tmp/Rtmp49kVVm/sphereL_158681.surf.gii /tmp/Rtmp49kVVm/sphereL_32492.surf.gii BARYCENTRIC /tmp/Rtmp49kVVm/to_read.surf.gii

ERROR: NAME OF FILE: to_resample.surf.gii
PATH TO FILE: /tmp/Rtmp49kVVm

and then it began those error messages about the zeros in the triangle array.

Also, I posted about this in the workbench github in case it is something on their end, but they indicate that they commonly work with very high res surfaces.

Best,
Steven

[damondpham]

Thanks again Steven for looking into this problem with us!

To recap the original problem: to resample the subject surface you've provided here, you will need to have a sphere at the same resolution (158681 vertices). But when I try to make such a sphere with the Workbench command "-surface-create-sphere", the resulting sphere has 158762 vertices, not 158681. So unless you already have a sphere with 158681 vertices, I'm not sure how you can resample this file.

I looked more into this follow-up problem, where it seemed like the Connectome Workbench cannot resample any large surface GIFTI file. Actually, it seems like it cannot resample any large surface GIFTI file in which the faces data is encoded with ASCII. ciftiTools::resample_surf uses gifti::write_gifti which, at the moment, can only encode integer data with ASCII. So a large GIFTI file can be resampled with the Workbench successfully if the faces data is encoded with e.g. GZipBase64Binary. But if it's read in with gifti::read_gifti and then written back out with gifti::write_gifti, like in ciftiTools::resample_surf, the new file will have faces data encoded with ASCII, and it will not be compatible with the Workbench. I've attached the document I used to show this. test_resampling_large.zip

I will look for a solution to the follow-up problem. Meanwhile, as the developers of the Workbench recommended, you could use -surface-resample more directly if you do have a sphere compatible with your high-res subject surface:

# Make the spheres you need
ciftiTools:::write_spheres("sphereL.32k.surf.gii", "sphereR.32k.surf.gii", 32000)
ciftiTools:::write_spheres("sphereL.164k.surf.gii", "sphereR.164k.surf.gii", 164000)
# Function to do resampling without reading & writing the file again, as in `resample_surf`
do_resample <- function(x, x_sph, y_sph, y_fname=NULL){
  yy <- tempfile()
  run_wb_cmd(paste(
    "-surface-resample", x, x_sph, y_sph, "BARYCENTRIC", yy
  ))

  # Return a `"surf"` object, or the name of the written file.
  if (is.null(y_fname)) {
    return(read_surf(yy))
  } else {
    file.copy(yy, y_fname)
    return(y_fname)
  }
}
do_resample(my_surf, sphere_originalRes, sphere_newRes)

[smeisler]

Great, thanks for the thorough investigation!

[coalsont]

Finally checked back on this. wb_command can operate on ascii-encoded large gifti files it generates:

wb_command -surface-create-sphere 180000 test_sphere_180k.surf.gii
wb_command -gifti-convert ASCII test_sphere_180k.surf.gii test_sphere_180k_ascii.surf.gii
wb_command -file-information test_sphere_180k_ascii.surf.gii

Perhaps something is going wrong in the library you are writing with, and your file ends up corrupted. It looks like we don't check for an insufficient number of items in ASCII encoding in workbench, so it may fill any missing parts of the file with 0 before it gets to the sanity check that you hit. Do a line count between the start and end tags of the triangles data array, and see if it makes sense, like so:

$ grep -n '<DataArray' test_sphere_180k_ascii.surf.gii
35:<DataArray Intent="NIFTI_INTENT_POINTSET"
179634:<DataArray Intent="NIFTI_INTENT_TRIANGLE"
$ grep -n '</DataArray' test_sphere_180k_ascii.surf.gii
179633:</DataArray>
538772:</DataArray>

538772 - 179634 = 359138, and this triangle array's Dim0 is 359120 (for 179562 actual vertices, since -surface-create-sphere just gets as close as it can), with the attributes, metadata, and formatting making up the extra 18 lines.

[damondpham]

Thank you Tim! I will do the line count check, and look into ciftiTools and gifti to see if we can fix this problem.

[damondpham]

Hi @coalsont,

Sorry for the delay with my response.

To restate the problem: If I write a large ASCII-encoded surface geometry GIFTI file with gifti::write_gifti, the resulting file cannot be resampled with -surface-resample. -file-information also fails. However, I can still read it in with gifti::read_gifti.

Here's what I see in the command line for a large non-ASCII-encoded surface geometry GIFTI file ("bigInflated.L.surf.gii") versus its ASCII-encoded counterpart ("bigInflated2.L.surf.gii"):

Non-ASCII large GIFTI file

PS the_dir> grep -n '<DataArray' bigInflated.L.surf.gii
48:<DataArray Intent="NIFTI_INTENT_POINTSET"
88:<DataArray Intent="NIFTI_INTENT_TRIANGLE"

PS the_dir> grep -n '</DataArray' bigInflated.L.surf.gii
87:</DataArray>
105:</DataArray>

ASCII large GIFTI file

PS the_dir> grep -n '<DataArray' bigInflated2.L.surf.gii
13:  <DataArray Intent="NIFTI_INTENT_POINTSET" DataType="NIFTI_TYPE_FLOAT32" ArrayIndexingOrder="RowMajorOrder" Dimensionality="2" Dim0="179562" Dim1="3" Encoding="GZipBase64Binary" Endian="LittleEndian" ExternalFileName="" ExternalFileOffset="0">
30:  <DataArray Intent="NIFTI_INTENT_TRIANGLE" DataType="NIFTI_TYPE_INT32" ArrayIndexingOrder="RowMajorOrder" Dimensionality="2" Dim0="359120" Dim1="3" Encoding="ASCII" Endian="LittleEndian" ExternalFileName="" ExternalFileOffset="0">

PS the_dir> grep -n '</DataArray' bigInflated2.L.surf.gii
29:  </DataArray>
1077401:  </DataArray>

I'm not sure how to interpret this result.

Here is the problematic GIFTI file: "bigInflated2.L.surf.gii". https://drive.google.com/file/d/1UwRAsGzG5G8eiMxQRJvqegyrto1noNfE/view?usp=sharing . Is truncation happening, or do you have any other ideas what the problem is?

Thanks for your help!

Damon

[coalsont]

Well, it wasn't truncated, but I found the issue in your file:

...
99999
99999
99869
1e+05
99869
99870
1e+05
1e+05
99870
100001
...

When reading an integer from text, we don't support scientific notation. I was worried that this was also causing a rounding error, but it looks like it only abbreviated exactly 100,000 this way (and presumably would do so for 200k, and 300k, and...).

[damondpham]

Thank you @coalsont !! I added a patch to gifti that seems to solve this problem by avoiding scientific notation.

[coalsont]

To be clear, I believe we do support scientific notation when the array is of type float32, it is just unexpected for integer types (because they need to remain exact).

[smeisler]

Hi, I am still getting this error using the most recent version of ciftitools.

[damondpham]

Yes, unfortunately gifti has not been updated CRAN, so the old version of gifti is causing the problem.

Could you try using the new gifti version?

devtools::install_github("muschellij2/gifti")

[smeisler]

The error persists on this version of gifti too.

[smeisler]

The error is different than the original error (Vertex used twice in one triangle).

Command

x <- resample_surf('/path/to/...._hemi-L_midthickness.surf.gii',10000,hemisphere=c("left"))

Error

While running:
/rdma/vast-rdma/vast/gablab/smeisler/workbench/bin_rh_linux64/../exe_rh_linux64/wb_command -surface-resample /tmp/RtmpHRZNCO/to_resample.surf.gii /tmp/RtmpHRZNCO/sphereL_227575.surf.gii /tmp/RtmpHRZNCO/sphereL_10000.surf.gii BARYCENTRIC /tmp/RtmpHRZNCO/to_read.surf.gii

ERROR: input surface has different number of nodes than input sphere

read_surf information of surface I am trying to resample

Vertices: 227575
Faces: 455146
Hemisphere: left

[coalsont]

Here is where that error text is, it is a simple sanity check that the input files look compatible:

https://github.com/Washington-University/workbench/blob/master/src/Algorithms/AlgorithmSurfaceResample.cxx#L150

Excerpt from the help info:

   wb_command -surface-resample
      <surface-in> - the surface file to resample
      <current-sphere> - a sphere surface with the mesh that the input surface
         is currently on
...

Check the contents of to_resample.surf.gii and compare to sphereL_227575.surf.gii. If those don't have the same number of vertices, wb_command -surface-resample can't deal with this, because the inputs are incompatible (having the same mesh is somewhat more specific than matching the number of vertices, but it is a start). Maybe the resample_surf function expects only freesurfer native mesh surfaces as inputs, and has no input option to use a different starting sphere? In most cases, you should only be resampling from the original freesurfer mesh version, to avoid multi-resampling loss.

[damondpham]

Thank you @smeisler for letting me know about this persisting problem, and thanks very much @coalsont for the guidance! I will look into this problem. I may have some follow-up questions, but in the meantime, @smeisler could you use the wb_command -surface-resample directly, with the original freesurfer spheres in the same resolution as the surface you are trying to resample?

If you don't have the compatible freesurfer spheres, I'm think that resampling will not be possible with the Connectome Workbench, nor will it be possible with ciftiTools.

[damondpham]

I re-read the thread above. To try to clarify: it seems like there were originally two problems.

1. gifti was writing large surfaces incorrectly.
2. ciftiTools was unable to create spheres that matched the resolution of the surface to be resampled.

I believe the first problem has been fixed. But the second problem persists, and cannot be fixed. I think the only solution is to use wb_command -surface-resample directly with the original freesurfer sphere in the matching resolution. Do you have that?

[smeisler]

Yes, I have the original FS spheres (/fsdir/sub-xxx/surf/lh.sphere, for example, or should it be lh.sphere.reg?). The syntax for -surface-resample is:

wb_command -surface-resample
      <surface-in> - the surface file to resample
      <current-sphere> - a sphere surface with the mesh that the input surface
         is currently on
      <new-sphere> - a sphere surface that is in register with <current-sphere>
         and has the desired output mesh
      <method> - the method name
      <surface-out> - output - the output surface file

What would I put in for <new-sphere>? (edit, I see now that I should be finding sphereL_10000.surf.gii. Trying to see how to get that since /tmp files are erased)

[smeisler]

Weird, I try the original x <- resample_surf('/path/to/...._hemi-L_midthickness.surf.gii',10000,hemisphere=c("left")) on a different subject, and I get the following error message instead:

Error in read_xml.character(file): Extra content at the end of the document [5]
Traceback:

1. resample_surf("/om4/group/gablab/data/inside/data//derivatives/fmriprep_22.0.0RC0/sub-inside7004/anat/sub-inside7004_hemi-L_midthickness.surf.gii", 
 .     10000, hemisphere = c("left"))
2. resample_gifti(gii_pre, gii_post, hemisphere = hemisphere, original_res = original_res, 
 .     resamp_res = resamp_res)
3. make_surf(original_fname, hemisphere)
4. read_surf(surf = surf, expected_hemisphere = expected_hemisphere, 
 .     resamp_res = resamp_res)
5. readgii(surf)
6. read_xml(file)
7. read_xml.character(file)

So it looks like this subject got further than the other, but is running into a later error.

[damondpham]

Hi @smeisler to create <new-sphere> you can use Connectome Workbench commands to create a sphere in the desired resolution. Or, you can an internal ciftiTools function ciftiTools:::write_spheres which wraps around the relevant Workbench commands. The source code is copied below for your reference:

write_spheres <- function(
  sphereL_fname, sphereR_fname, resamp_res, 
  write_dir=NULL) {

  sphereL_fname <- format_path(sphereL_fname, write_dir, mode=2)
  sphereR_fname <- format_path(sphereR_fname, write_dir, mode=2)

  resamp_res <- format(resamp_res, scientific=FALSE)

  # Make left
  run_wb_cmd(
    paste("-surface-create-sphere", resamp_res[1], sys_path(sphereL_fname)), 
  )

  # Make right
  if (length(resamp_res)==1 || length(unique(resamp_res))==1) {
    run_wb_cmd(
      paste("-surface-flip-lr", sys_path(sphereL_fname), sys_path(sphereR_fname)), 
    )
  } else {
    sphereR_temp <- paste0(tempfile(), ".surf.gii")
    run_wb_cmd(
      paste("-surface-create-sphere", resamp_res[2], sphereR_temp), 
    )
    run_wb_cmd(
      paste("-surface-flip-lr", sphereR_temp, sys_path(sphereR_fname)), 
    )
  }

  # Set structure
  run_wb_cmd(
    paste("-set-structure", sys_path(sphereL_fname), "CORTEX_LEFT"), 
  )
  run_wb_cmd(
    paste("-set-structure", sys_path(sphereR_fname), "CORTEX_RIGHT"), 
  )

  invisible(list(sphereL_fname=sphereL_fname, sphereR_fname=sphereR_fname))
}

With the proper <current-sphere> and <new-sphere>, you should be able to resample your file!

[smeisler]

I was able to resurface the sphere and surface using the wb-commands directly. Minimal example:

ciftiTools:::write_spheres('/path/to/freesurfer/$SUB/surf/lh.sphere',
                           '/path/to/freesurfer/$SUB/surf/rh.sphere',
                           10000, write_dir='resampled')

mris_convert /path/to/freesurfer/$SUB/surf/lh.sphere mris_convert /path/to/freesurfer/$SUB/surf/lh.sphere.gii

(REPEAT WITH RH, INPUT RESAMPLED SPHERES)

wb_command -surface-resample /path/to/fmriprep/$SUB/anat/$SUB_hemi-L_midthickness.surf.gii /path/to/freesurfer/$SUB/surf/lh.sphere.gii /path/to/resampled/lh.sphere BARYCENTRIC /outdir/lh_surf_resampled.surf.gii
(REPEAT WITH RH)

Will I have to resample my func.gii before inputting to BayesGLM if I want to completely bypass the built in resampling method?

[damondpham]

That's great you were able to! Unfortunately yes, you would need to resample separately as you did here, and then call BayesGLM after. But I will look into implementing this type of resampling within BayesGLM.

[coalsont]

Yes, I have the original FS spheres (/fsdir/sub-xxx/surf/lh.sphere, for example, or should it be lh.sphere.reg?). The syntax for -surface-resample is:

wb_command -surface-resample
      <surface-in> - the surface file to resample
      <current-sphere> - a sphere surface with the mesh that the input surface
         is currently on
      <new-sphere> - a sphere surface that is in register with <current-sphere>
         and has the desired output mesh
      <method> - the method name
      <surface-out> - output - the output surface file

Sorry I didn't see this earlier. The current-sphere and new-sphere need to be "in register", so if one of them is a native mesh sphere and the other is a standard template sphere, one of them will need to be the "after registration" sphere, otherwise your data won't line up with other subjects using that standard mesh. However, if you are resampling to a standard mesh other than fsaverage, then freesurfer's sphere.reg won't give you a correct result, either (I think you were aiming for fs_LR 32k originally, which is not fsaverage). For this case, the HCP pipelines have fs_LR-deformed_to-fsaverage spheres in the global/templates/standard_mesh_atlases/resample_fsaverage/ folder - using those with the subject's sphere.reg should give you correct results on the fs_LR mesh. For a more complete procedure, including using the adaptive barycentric method to avoid SNR penalties when downsampling (though this may not matter for anatomical measures), see HCP users FAQ #9, and look at section B in the pdf.

Edit: I will also note that the HCP pipelines use MSMSulc and MSMAll registration rather than freesurfer registration, which result in less mesh distortion, and particularly for MSMAll, better functional area alignment across subjects as measured by task data. Running MSMSulc may be fairly straightforward on data that hasn't been processed with the HCP pipelines, but running MSMAll would be more of a challenge (it needs a fair number of files, and expects them in the HCP pipelines organization).