Imperfect overlap of WM and GM mask, especially in in the lower thoracic and lumbar levels

[jcohenadad]

@Nilser3 noticed some issues with the GM mask PAM50_gm.nii.gz in version 6.5:

Issue 1: Mismatch between GM mask and T2* template

Issue 2: Slight mis-registration between the T2*w and the T2w template

Issue 3: Some pixels far away from the GM mask have non-zero values

[Nilser3]

Issue 4: Mismatch of the GM / WM masks with PAM50_t2s template

EDIT by @joshuacwnewton: Since the WM mask is generated by subtracting the GM mask from the SC mask), then any GM problems will persist in the WM mask. In other words, I think the problems with the WM mask come directly from the problems with the GM mask.

This means that this sub-issue is a duplicate of Julien's report of "Issue 1" (problems with the GM mask). To make this issue easier to follow, I'm collapsing it as "Duplicate".

Details

Looking in more detail between PAM50_gm.nii.gz , PAM50_wm.nii.gz and PAM50_t2s.nii.gz, we note that the mismatch also persists in all the PAM50_t2s template (ie. Level C2):

PAM50_gm.nii.gz and PAM50_t2s.nii.gz

PAM50_wm.nii.gz and PAM50_t2s.nii.gz

• Version of PAM50 : SCTv.7.0

To see the mismatch (in the whole PAM50_t2s.nii.gz ), we have to adjust the contrast.

Related issues

#33

Tagging @valosekj

[Nilser3]

Comparison with May 2016 AMU15 template files

Checking if the mismatch presented above comes from previous versions, it was verified that the AMU15 template (files from May 2016) did not present this mismatch (same axial slice = 933):

Large screenshots

Template and GM mask

Template and WM mask

Template and GM mask in lower levels

Template and WM mask in lower levels

I imagine that in the conversion from AMU15 to PAM50 there was some post-processing steps that changed the GM and WM mask?

Tagging @vcallot

[vcallot]

@Nilser3 , thanks for illustrating the issues.
@jcohenadad , @valosekj : to complete what Nilser has indicated : the last AMU15 version we (Marseille) provided was for the generation of the MNI-POLY-AMU version3 (release in april 2016 by Benjamin), illustrated above.

[jcohenadad]

thank you @vcallot and @Nilser3 for the additional context-- we will discuss it in our next SCT meeting and come up with a plan

[joshuacwnewton]

I imagine that in the conversion from AMU15 to PAM50 there was some post-processing steps that changed the GM and WM mask?

Note: The CHANGES.md file contains some commands that seem to have been used to register AMU -> PAM50 in July 2016.

Using this as a starting point, I traced the history of the PAM50 template:

Release	Image	Notes
r20160716: First release within the PAM50 GitHub repo. Associated with "# register AMU->PAM50" code.		Segmentation quality is very poor in lower levels. But, the issues are different than what is shown in this issue.
r20160906: Missing CHANGES.md entry. Linked issue mentioned in the associated commit (43813bf - #968) suggests C1 extension.		Improved segmentation quality. GM+WM now sum to 1. First release with the "shifted lower" mis-registration issue.
r20181214: From CHANGES.md: "Aligned the template with the ICBM152"		Scanner coordinate update. No changes to GM/WM segmentation shape. However, GM+WM no longer sum to 1.0 (!!!).
r20231024: From CHANGES.md: "Add script that symmetrizes images."		Scanner coordinate update. No changes to GM/WM segmentation shape.
r20250422: Latest release.		Scanner coordinate update. No changes to GM/WM segmentation shape. GM+WM sum to 1 again (!!!).

So, I think the key takeaways are this:

• The shape of the GM/WM segmentation has not changed since September 2016.
• The July 2016 sct_ commands seem to be irrelevant, and instead the issues come from whatever happened afterwards between July and September 2016. Unfortunately, however, the PAM50 changelog does not include an entry for the September release to explain what commands were run.
• The commit that updated the PAM50 URL in SCT is spinalcordtoolbox/spinalcordtoolbox@1284a9b. This commit doesn't include any code either.
• There are some DEV: commits a few weeks prior that update a script called create_atlas.m. Diffing create_atlas.m between July 11 and Sept 1, significant changes were made in these DEV: commits: https://www.diffchecker.com/CkBZxGqZ/.
• I think especially of note is commit spinalcordtoolbox/spinalcordtoolbox@dc02834:

  "now generates WM, GM and cord (binary) global atlases that are the exact sum of WM/GM individual atlases"

  • This lines up with the fact that the Sept 2016 release changes the shape of the GM+WM masks, and that the masks summed to 1 starting with that release as well.
• So, it seems create_atlas.m is the source of the issues found in the files within our current PAM50 template.

---

As far as actually fixing this issue, I would guess that we would want to recreate the same environment as when the create_atlas.m script was run in September 2016, then tweak the parameters. However, this might be tricky:

• The create_atlas.m script references data on Julien's local drive. It's not clear to me where this data may have been stored long-term.
• The script relies on specific MATLAB libraries as well as c3d and ANTs as dependencies, and so we'd need to be careful to find the versions used in 2016 if we wish to recreate the original conditions (which we can then tweak). If we use more modern versions of these dependencies, there is a chance that upstream code changes could disrupt the results of the template generation (breaking API changes, numerical differences, etc.).
• The fact that the GM/WM segs are so closely tied to the GM/WM atlas suggests that if we change the segmentations, we'd also have to change the atlas, and that seems like a much more consequential change (given our atlas-based analysis).

I think that's as much as I can do on my own without further input, so if there's something I've missed (or if there are ideas on how to move forward), please feel free to chime in. 😄

[jcohenadad]

The create_atlas.m script references data on Julien's local drive. It's not clear to me where this data may have been stored long-term.

@joshuacwnewton unfortunately i cannot see the data on my own servers (because of https://github.com/neuropoly/computers/issues/683#issuecomment-3014221569), but i think there is a dataset called PAM50 or something... maybe this is now called whole-spine if there is no other dataset

The script relies on specific MATLAB libraries as well as c3d and ANTs as dependencies, and so we'd need to be careful to find the versions used in 2016 if we wish to recreate the original conditions (which we can then tweak). If we use more modern versions of these dependencies, there is a chance that upstream code changes could disrupt the results of the template generation (breaking API changes, numerical differences, etc.).

i think we can go with more recent versions, and make sure results "look good". If we have minor numerical changes, this is not the end of the world.

The fact that the GM/WM segs are so closely tied to the GM/WM atlas suggests that if we change the segmentations, we'd also have to change the atlas, and that seems like a much more consequential change (given our atlas-based analysis).

the atlas would need to be changed, yes, because it is constructed by create_atlas.m. Hopefully the change won't be too drastic-- i expect minor numerical changes

[joshuacwnewton]

Thank you for the details! OK, I can take a stab at finding the data + writing a set of steps needed to recreate the environment. :)

[joshuacwnewton]

I've been able to successfully run create_atlas.m!

To demonstrate, I created a branch jn/32-regenerate-atlas-and-gm-masks that contains both the necessary data files, and an updated copy of create_atlas.m that fully explains how to set up the necessary dependencies (as well as the versions I tested):

PAM50/scripts/create_atlas/create_atlas.m

Lines 10 to 37 in 61dfb41

	%----------------------------- Dependencies -------------------------------
	% 1. MATLAB: Ensure both Image Processing Toolbox and Signal
	% Processing Toolbox are checked during installation.
	% (R2025a tested)
	% See: https://www.mathworks.com/help/install/ug/install-products-with-internet-connection.html
	% 2. FSL: Download and install FSL (FMRIB Software Library).
	% (v6.0.7.17 tested)
	% See: https://fsl.fmrib.ox.ac.uk/fsl/docs/#/install/index
	% 3. c3d: Download and install c3d (Convert3D).
	% (Newest version failed (v1.4.2) -> Older version succeeded (v1.0.0))
	% See: https://sourceforge.net/projects/c3d/
	% 4. ANTs: Download and install SCT, which contains the necessary ANTs scripts.
	% (SCT v7.0 tested)
	%
	% - Additionally, ensure the following scripts are in the working directory:
	% * m_normalize.m : normalization function
	% * dnsamplelin.m : function for downsampling by computing mean value for
	% each region
	% * m_linear_interp.m
	%
	% - Additioanlly, ensure the following data files are present in ../raw_data:
	% * PAM50_wm.nii.gz
	% * atlas_grays_cerv_sym_correc_r6.png
	% * atlas_grays_cerv_sym_correc_r6_label.txt
	% * mask_grays_cerv_sym_correc_r5.png
	%
	% You can get these files from the following link (preserved by GitHub):
	% https://github.com/spinalcordtoolbox/spinalcordtoolbox/tree/1284a9b09766916cd685c77bfdc3bb2f68c07a98/dev/atlas/raw_data

The output files have slightly different "stray nonzero voxels" than what I observed in our current template, but the general shape is basically the same as our current template:

There are also some issues with the immediate output:

• GM+WM don't sum to 1 (as observed in previous copies of the template)
• Header is not correct (large offset between the generated files and the PAM50 T2* template image.

Still though, I think we have a good starting point here for tweaking the results. 😃

Note

The current create_atlas.m procedure takes quite a while (20-30m on my machine), so we may want to be smart with how we tweak the parameters given the time investment involved in each attempt to regenerate the atlas. I may need to return to this issue once v7.1 has been completed.

[jcohenadad]

Amazing! I'm glad you were able to reproduce the results Joshua. I'm not soooo surprised, though, that's one of the (rare) good thing about Matlab: reproducibility and cross-platform compatibility.

The current create_atlas.m procedure takes quite a while (20-30m on my machine), so we may want to be smart with how we tweak the parameters given the time investment involved in each attempt to regenerate the atlas. I may need to return to this issue once v7.1 has been completed

I wouldn't worry about that too much, given how unfrequent we generate this atlas-- even for prototyping. Unless there is an obvious and easy way to accelerate things of course 😊