Add older processing scripts to clarify provenance of PAM50 atlas files #47
Conversation
This should be everything that's needed to run the script. :)
Previously it was too difficult to tell what part of the script you were in (since it takes a long time, it's important to know what steps we're on).
These files probably shouldn't be committed to the PAM50 repo, but I'm adding them to this branch just to make reproduction easier for anyone else checking out my branch.
Making it more discoverable -- though this doesn't fully document the exact procedures used. Instead, the CHANGES.md file should be used for this purpose.
This reverts commit 4b85d81.
There was a problem hiding this comment.
I tried tracing the script, setting breakpoints, examining the internal variables during this script... and there were a few parts that stood out to me?
But, ultimately, debugging MATLAB scripts is tricky! I don't quite understand the save_avw function (and how to save arrays as images). It's also quite time-consuming to get to the point where these steps would be run. So, I don't feel confident I've identified the right bits of code. 😓
Still, hopefully this is a useful starting point @valosekj?
| % Apply tranform to the PVE tract files | ||
| for label = length([label_left, label_right])+1:length(label_values) | ||
| tract_atlas = [ 'tract_atlas_' num2str(label)]; | ||
|
|
||
| cmd = ['isct_antsApplyTransforms -d 2 -i ' tract_atlas ext ' -o ' tract_atlas suffix_ants ext ' -t ' Warp_atlas ' ' affine_atlas ' -r ' templateci_slice_ref_thresh ext]; | ||
| disp(cmd); [status,result] = unix(cmd); if(status), error(result); end, disp(result) | ||
|
|
||
| tract_reg_g = [ 'tract_atlas_' num2str(label) suffix_ants]; | ||
| temp_g = read_avw(tract_reg_g); | ||
|
|
||
| % Replace isolated values with the mean of the adjacent values | ||
| for i = 2:size(temp_g,1)-1 | ||
| for j = 2:size(temp_g,2)-1 | ||
| test = (temp_g(i,j)==temp_g(i-1,j)) || (temp_g(i,j)==temp_g(i,j-1)) || (temp_g(i,j)==temp_g(i+1,j)) || (temp_g(i,j)==temp_g(i,j+1)); | ||
| if(~test) | ||
| temp_g(i,j) = (temp_g(i-1,j) + temp_g(i+1,j) + temp_g(i,j+1) + temp_g(i,j-1))/4; | ||
| end | ||
| end | ||
| end | ||
|
|
||
| tractsHR{label}(:,:,num_slice_ref) = temp_g; | ||
| end |
There was a problem hiding this comment.
These lines here are the first mention of PVE in the create_atlas.m script.
When I look at the files we're working on (e.g. tract_atlas_XY_reg.nii.gz), we can see that at this point, we're still working with the digitized, high-resolution Grey's Anatomy masks registered to the PAM50 space, without any shape transformations:
The code above:
- Iterates through voxels in each of the 2D tract images
- Computes
test(which checks if any of the adjacent voxels are equal to the current voxel) - If no adjacent voxels are equal, then a simple mean is computed. In practice, this means that the mean is only computed for the edges of each tract mask.
- This is kind of similar to a 3x3 Gaussian kernel in that it performs a weighted average of the neighbourhood, but the weights are strange? Instead of something like this:
What's actually used is this:
Additionally, these values are not immediately written to disk, but instead stored in-memory in a variable called tractsHR. I presume this means "high resolution".
Ultimately, I suspect that this averaging operation is a red herring, and doesn't have a significant effect on the PVE.
| %% Downsampling and partial volume computation | ||
| disp('*** Downsampling computed slices for each label... ***') | ||
| max_indx = max(z_disks_mid(:)); | ||
| for label = 1:length(label_values) | ||
| disp(['LABEL #: ', num2str(label), '/', num2str(length(label_values))]) | ||
| for zslice = 0:max_indx | ||
| numSlice = zslice+1; | ||
| tracts{label}(:,:,numSlice) = dnsamplelin(tractsHR{label}(:,:,numSlice),interp_factor); | ||
| end | ||
| end |
There was a problem hiding this comment.
I am wondering if this downsampling is what actually introduces the "soft" effect? 🤔
I'm in the process of recreating the T2*w template, and then the WM/GM atlas, so I will have to revamp that entire script (which I could port to Python while doing that), so what I'm trying to say is that in the end I could take care of this |
|
Thanks for digging into it, @joshuacwnewton! |
Gotcha-- AFAIR, no Gaussian filter was applied. The atlas was "simply" built by partitioning high res spinal cord mask into smaller non-overlapping regions, then downsampled, resulting in PVE information (no Gaussian smoothing). Then, slight non-linear deformation slicewise. |
|
That lines up with my understanding too after reading through the script. Glad to hear it. :) |
|
I see! Thanks for the details. |
3 participants
Description
Prior to my investigation in #32 (comment), it was unclear which scripts were used to generate the WM atlas files.
This PR adds what I believe to be the necessary files (though there are multiple "atlas creation scripts" so I may have the wrong one). This PR also adds links to the README pointing to the archival locations of the scripts.
We can use this as a baseline for further re-generation of the atlas (e.g. for #41).
Discussion
We can probably safely drop commit 4b85d81 from this PR, since we don't really want to commit unnecessary binary files to the repo, and the files are also backed up in the
spinalcordtoolboxrepo (and linked to in the "Depencencies" section of thecreate_atlas.mscript).But, I wanted to add them initially just for ease of access...
Related Issues/PRs
Fixes #46.
Part of #32.