Question on output of metric-gradient

[Jason Davis]

What should the size/output of the metric gradient file? For instance, if I have a similarity matrix of roughly 10K vertices (10k X 10k), my gradient files, produced by -metric gradient in the workbench are also the same matrix (10k X 10k)? Shouldn't the gradient values be a 10k vector? If not, do I need to average within the matrix in some way? How can I interepret the gradient of the matrix? Please help, as I am trying to recreate a version of the Gordon parcellation and am a little stuck on this step.

Jason

[Tim Coalson]

-metric-gradient uses a surface file and outputs the spatial gradient magnitude of each map independently.  The optional vectors output gives 3 maps per input map, which are the x, y, and z components of the gradient vectors (while the gradient vectors are always in the tangent plane, the plane normal is different at each vertex, so 3D is used for consistency of representation).

Yes, the -metric-gradient output with a 10k x 10k input will also be a 10k x 10k file, because each input map is treated independently.  If you want a command that averages the gradients across maps after computing them, look at -cifti-correlation-gradient.  We do not have a gradient-like command that treats the input maps as non-independent.  Also note the -double-correlation option, and that Gordon used 32k.  10k doesn't represent the geometry very well in adult humans, so the surface gradient will be less accurate than 32k.  However, Gordon didn't use wb_command's surface-based gradients for this step, and involved a fair amount of custom matlab code, including the watershed step.

Tim

[Jason Davis]

Hi Tim,

Thanks for responding so quickly. I appreciate your answer. So then if I were using GIFTI files, would a similar implementation be to get the gradient maps and then average across columns to get a value for each row (leading to a 10k by 1 column matrix)? If the rows represent physical locations, then if I understand correctly, couldn’t I average each gradient maps (each column) to get an average?

I also saw previous iterations in Wig (2014) use the -metric-gradient-all function in Caret 5.65. I wanted to take a look at that function but I couldn’t load the software. Is there any dramatic differences between HCP’s metric-gradient and this other Caret function?

Based on my goal (of trying to recreate the Gordon parcellation scheme as closely as possible), what seems like the most reasonable approach to match their work? I can’t really create CIFTI files since I don’t have any subcortical data.

The -cifti-correlation-gradient function asks one to take the gradients to each of the rows and then averages them. The -cifti-gradient function also seems to average the output of the gradient maps as well, assuming that you input this as a parameter. Would my solution of averaging the maps after creation suffice potentially?

Best,

Jason

[Jason Davis]

Hi Tim,

It also seems like the correlation values were done on each surface vertex compared with every other vertex across hemispheres, and then the similarity values were done by extracting the similarity values across all of the hemispheres and then producing a specific hemisphere matrix. 

Were the initial correlations done for each hemisphere separately (i.e. 32k X 32k correlation matrices for each hemisphere) or inclusive of data for both hemispheres (64k X 64k matrix)? 

If they included data from both hemispheres, how can there be hemisphere specific similarity matrices developed? Did they just extract the quadrant of the matrix that corresponds to each hemisphere? 

I am having trouble understanding the details of the pipeline. I don’t know if you can help explain it to me, but if you could, that would be greatly appreciated. Thank you!

Best,

Jason

[Tim Coalson]

The cifti format doesn't require you to include subcortical data, but our pipeline outputs would include the subcortical data (because it expects the volume files to have acquired the entire brain).  However you do it, you should exclude the medial wall from the cortical processing (with an roi to -metric-gradient or using a medial-wall-excluded cifti space).

Yes, you can do the steps separately and get the same answer (with more intermediate files and disk IO).  The caret5 metric gradient algorithm is very similar, though the workbench version is somewhat better at not being as biased by curvature.

The Gordon paper seems to say they did each hemisphere by itself for both correlation steps.  We would use the 91k grayordinates (both hems and subcortical) for the first correlation, and then the single-structure (e.g., single-hemisphere) correlation between the 91k-long maps for the second correlation step - this is what -cifti-correlation-gradient -double-correlation does.

Tim

[Jason Davis]

Hi Tim,

Thank you for your help. I am still a bit confused on the -roi parameter to specify in the -metric-gradient command? What do I have to specify in the actual metric to exclude the medial wall? We have a masking function that excludes vertices from the medial wall for each hemisphere? Would I have to create a mask file that maybe denotes where to take the gradient of (as in 1s for non-medial wall and 0s for medial wall)? It says it is in terms of area but could you do it by vertex number?

[Glasser, Matthew]

Yes, that is the idea.  1s in the non-medial wall cortex and 0s in the medial wall.

Matt.

--
You received this message because you are subscribed to the Google Groups "HCP-Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hcp-users+...@humanconnectome.org.
To view this discussion visit https://groups.google.com/a/humanconnectome.org/d/msgid/hcp-users/9fea902c-5177-41ef-a8b8-68cef89d26f6n%40humanconnectome.org.

 

---

The materials in this message are private and may contain Protected Healthcare Information or other information of a sensitive nature. If you are not the intended recipient, be advised that any unauthorized use, disclosure, copying or the taking of any action in reliance on the contents of this information is strictly prohibited. If you have received this email in error, please immediately notify the sender via telephone or return mail.

[Tim Coalson]

In the "-roi" option's description, the word "area" is being used in the generic sense of "arbitrary bounded region of the surface".  It does not refer to surface area, nor specifically to functional/architectonic areas.  It just takes a .func.gii file with values of 1s and 0s.

Tim

[Jason Davis]

How can I adjust the surface file to make this computation possible? I tried using the -roi parameter and the masked ROI in this case (LH - 9675, RH - 9666 vs. unmasked - 10242) will not match the number of vertices in the original surface file (10242). I tried masking the surface file for the first data array but I am not sure what to do with the second data array in the surface file?

[Tim Coalson]

Don't create a shorter array, every vertex needs to be represented in all gifti files.  The roi file needs to be 10242 vertices, just like your surfaces and data.  You shouldn't be trying to modify any .surf.gii files, that would require a lot of extra complication.

If you aren't using fs_LR template conventions, then the only easy option is to resample the subject's freesurfer medial wall to 10242 (and logically negate if necessary so that cortex is where the 1s are).  If you are using fs_LR conventions, you could instead resample our standard cifti medial wall files to 10242, which are the .shape.gii files here:

https://github.com/Washington-University/HCPpipelines/tree/master/global/templates/91282_Greyordinates

ROI files should generally be resampled with interpolation and then thresholded at 0.5 (assuming the ROI was 1s and 0s) for the best quality and to avoid accidental dilation-like effects.

Tim

[Jason Davis]

Hi Tim,

The challenge is that I am working with the onavg template, created by Feilong Ma in 2024 that more evenly samples the cortex. All of the files in the dataset have been resampled into that space already, so it would be difficult to go back and resample them back into fsLR space. The current data uses the ico32 resolution which contains 10242 unmasked vertices for each hemisphere (9675 - LH, 9666 - RH). 

I had already tried to modify the surf.gii files to only extract coordinates from the masked vertices in each hemisphere and then keep the triangles from vertices only included in the masked array, so that all of the triangles on the surface that are in the medial wall are not considered. Could you please explain why I shouldn’t do this? 

I had begun to run some code using these shorter files and it was able to work, although just because something works doesn’t mean it’s right, of course. 

I’m not really sure other ways to exclude the medial wall except for after the gradient calculation. I don’t think it makes sense to resample the freesurfer medial wall. It seems like our group is mainly focused on using the onavg template, which is what makes this more difficult. It doesn’t seem straightforward to fix the number of vertices accordingly. 

Any thoughts or suggestions? It seems like a real sticking point. 

Thanks for your help!

Jason

On Jun 3, 2026, at 8:53 PM, Tim Coalson <tim.c...@gmail.com> wrote:



[Tim Coalson]

I am not aware of onavg recommending that people delete vertices from their surface files.  Their python examples do appear to use fewer indices for their data than for their surfaces (when "mask=True", but that also isn't what we mean when we say "mask" - to us, "mask" means "set some data to zero, but don't alter the size of any arrays"), though it isn't clear where they actually store the mapping back to the full surface (which appears to be present in their example's plots, as the medial wall is explicitly drawn with a different color), so it appears that they do not delete the vertices from their surfaces.  Our tools don't use their "mask=True" convention, at least not in gifti format, and will not be able to display or operate on such incomplete arrays with the full surfaces.  The ability to save such collapsed data into gifti format is a fairly unfortunate hazard.  I would imagine that their python code has a function somewhere to turn the collapsed data back into the full 10242 vertices it logically represents - but please don't put NaNs into the medial wall, just use 0s instead (NaN propagates through any floating point operation, so if a spatial operation doesn't use exactly the right ROI, a bunch of data near the medial wall gets replaced with NaN, which quickly invades the rest of the surface through later steps, even if those steps are fully correct).

The integers in the triangles array refer to the indices of the coordinates array.  Removing indices from the coordinates array without very carefully altering the triangle information will make your surface broken, and everything that uses spatial neighbor information will break and give the wrong answer (if it doesn't error first).  Even done correctly, altering your surfaces this way will likely make them incompatible with the python code onavg used, and your data will be superficially incompatible with any dataset that used the onavg spheres with preexisting tools outside of python.  Our surface-data resampling commands also expect the full sphere surface, and for the data to have the same number of vertices as the sphere.  Gifti files are usually internally compressed anyway, so any data file that has zeros in the medial wall may already be using little to no disk space to represent the data values of those vertices.

The cifti format, on the other hand, does in fact (usually) exclude medial wall vertices from its data files (but the .surf.gii files aren't altered, so the cifti header contains the mapping back to the full surfaces, and volume if applicable).  Consequently, the way that most spatial cifti processing works...is by expanding it back out to all vertices for one hemisphere at a time, doing the operation, and then recompressing the result, because anything else would be a lot more difficult for minimal efficiency gains.  The original intention of cifti was to enable correlation across all brain structures without including voxels and vertices that don't have usable data.

ROI options already exist for all spatial processing of gifti data in wb_command.  We even have an roi option for surface-based resampling, but like all the others, this only prevents nonsensical data from being used, it doesn't remove vertices from the output array.

In short, our tools are built around the assumption that the subject data on standard mesh is the same dimension as the standard spheres and other surfaces, without removing any vertices.  Using "mask=False" convention when saving onavg python data to gifti format will save you, and those working with your data in the future, a lot of time and effort.

Tim

[Jason Davis]

Hi Tim,

Thank you so much for your clarification. It makes sense that the tools rely on the subject data being the same dimensionality as the standard surfaces without removing any vertices. I believe the onavg version removes those vertices in the timecourse but the plotting function automatically knows to set those values to 0 if the data array input equals the masked number of vertices. So for example, if I had a timecourse of (1, 20484) which is 10242 vertices unmasked for each cortex and you input a data array of (1, 19341), then it recognizes the shape and plots the data as such.

Removing vertices from the output array doesn't seem like the best move. I see that Workbench has a -metric_mask function that can do this for me. I could then input the mask for the onavg space and the mask function will input 0s into the medial wall. As such, the data should be able to retain the same number of vertices and then I can use the -metric gradient function to calculate the gradient. That way, the medial wall is already being excluded from calculation but the similarity file input for the gradient retains the same number of vertices as the surface. Could this work as a possible solution?

If so, would it be best to mask the timecourses themselves using -metric_mask before doing the RSFC similarity calculation (as in the similarity calculation does not include the medial wall and this file continues as input into the gradient function)? Or should I conduct the RSFC calculation with the medial wall timecourses included and then mask the RSFC similarity matrix file? If the cifti files already exclude the medial wall, it seems like the best solution would be the former (masking the data before RSFC similarity calculation to ignore the medial wall as a factor for RSFC calculation?

This has been very helpful and clear. Thank you again for your help!

Jason

[Jason Davis]

Would I potentially need to mask both files using -metric-mask?

[Tim Coalson]

If you just -metric-mask, then -metric-gradient will see zeros in the medial wall, and will therefore compute a strong gradient exactly on the edge of the mask ROI.  The -roi option of -metric-gradient itself is what you need, to tell the gradient computation itself to ignore the data in those vertices, regardless of whether you replace the nonsensical data with zeros in advance.

Zeroing out timecourses will cause correlation to output NaN values, so that isn't a great idea.  Cifti format is how we usually do correlations, and all of our cifti commands automatically use the information of excluded vertices to restrict spatial operations from using the data that isn't in the cifti file.  Usually we use cifti to put both hemispheres into one cifti file, but if you are trying to exactly match a specific prior implementation that didn't correlate across hemispheres, it is possible to make single-hemisphere cifti files.  Look at wb_command -cifti-create-dense-timeseries and similar.

Tim

[Jason Davis]

Hi Tim,

I realized that you recommended two options: (a) implementing the -roi parameter to -metric-gradient or (b) using a medial-wall excluded CIFTI space. It seems like the -roi parameter option is quite cumbersome, so I am trying to create the respective CIFTI files. It seems as if the paper says that they correlated across hemispheres, although I could be mistaken. Would this mean that it is correlated to every other surface vertex and subcortical voxel (within the respective hemisphere) or is this meaning that it is correlated to every surface vertex and subcortical voxel (including those from other hemispheres). If they are producing these 32k maps, it seems like it was done based on the specific hemisphere rather than across hemispheres. This is the text I am referring to:

For each subject, the time course of each surface vertex was correlated with that from every other surface vertex and subcortical voxel in
the CIFTI space. Each correlation map was transformed using Fisher’s r-to-z transformation. For each hemisphere, the subject’s RSFC map similarity matrix was created by calculating the pairwise spatial correlations between all vertex’s RSFC correlation maps, producing a 32k × 32k matrix.

Additionally, can I use the -cifti-correlation-gradient or the -cifti-gradient function on single hemisphere cifti files? Or would it need to include both hemispheres and/or subcortical data? I don't believe we have subcortical data so this may be an issue. From your previous message, it seems like the pipeline expects subcortical data.

Jason

[Jason Davis]

I assume I can't use the -metric-gradient on the CIFTI file, correct? I assume you would have to use the CIFTI version of the function. I was wondering whether I can use the same function for each hemisphere, assuming I couldn't use the -cifti-gradient or -cifti-correlation-gradient functions accordingly (especially if it requires subcortical data).

[Tim Coalson]

Correct, cifti data requires -cifti-* commands.  All -cifti-* commands should accept single-hemisphere cifti files if you want to do it that way.

I'm not sure how I previously got the impression they hadn't correlated across hemispheres in the first step.  The "32k x 32k" part appears to refer to the second correlation step, not the first.  It sounds like all three steps are what the -cifti-correlation-gradient command with the -double-correlation and -fisher-z-first options does internally.  So, it sounds like you should create a combined cifti timeseries file (both hemispheres, and subcortical) if you want to reproduce what Gordon did.

I will mention again that 10k is rather low resolution, even for using fMRI data, such that the inaccuracy of the surface geometry will probably have an impact on the gradient calculation.  The cifti commands that use a lot of memory, including -cifti-correlation-gradient, have options to limit the memory usage (though some of them may be more approximate than others), which should make it possible to run at 40k on a modest computer, given enough extra time (correlation is symmetric, so computing the full matrix in memory can do just half the computation, but if you can only hold small bits of it at once, you have to compute nearly the whole thing - this gets even more complicated when there are two correlation steps, so try not to crank the limit too far down, or you will be waiting a long time).

Tim

[Jason Davis]

After computing the gradient using -cifti-correlation gradient, I do get the gradients, but it seems like they averaged across gradient maps, and -cifti-correlation-gradient doesn't seem to give you the actual matrix. It only gives you the final output. There's no option to get the gradient maps using this function before averaging, correct? I would have to use the -cifti-gradient to complete this.

Additionally, what is the most straightforward way to convert a RSFC similarity matrix into a CIFTI file? I calculated the original 10242x10242 similarity matrix and then input into a CIFTI using -create-dense-scalar with the ROI information, but it produces a 10242x9675 (masked LH) and 10242x9666 (masked RH) matrix, which is incorrect. If I input this into the -cifti-gradient function, will it be able to extract the ROI information and produce the correct information? I am not sure what is the easiest way to go about this.

Jason

[Tim Coalson]

That complication is a substantial part of why -cifti-correlation-gradient doesn't out the intermediate files.  For our 32k surface data, the first correlation is a 91k x 91k matrix, the second correlation and the gradient are pairs of 29k x 29k matrices of slightly different sizes (ignoring subcortical structures for the moment), but the final output is a single averaged map, which is much simpler to output.

Cifti format, unlike gifti, generally should have the medial wall excluded, because cifti format can store the information of what vertices were kept.  You don't want to use the gradient of maps of "what is similar to this medial wall vertex", so those maps should also be excluded.

You can do it the extra file IO and large files way with separate commands, but it is more complicated and slower, and will produce the same final map when done right.  Why do you want the large intermediate matrices?  I would probably start with both-hemisphere cifti timeseries with the medial wall excluded, -cifti-correlate, abuse -cifti-resample to separate out the single-hemisphere rectangles of the both-hemisphere matrix directly to cifti format, and then -cifti-correlate, -cifti-gradient, and -cifti-reduce.

Tim