Add older processing scripts to clarify provenance of PAM50 atlas files

README.md

@@ -4,3 +4,8 @@ PAM50 template and atlas of the white and gray matter spinal cord. For more deta
 
 ![Screen Shot 2021-08-20 at 12 01 54 PM](https://user-images.githubusercontent.com/2482071/130217021-d0064253-a540-4b39-b426-d1fe731ceb4c.png)
 
+## Useful links 
+
+- [`PAM50/CHANGES.md`](https://github.com/spinalcordtoolbox/PAM50/blob/master/CHANGES.md): Changelog of the PAM50 template and atlas files. Includes details about how certain files were generated.
+- [`PAM50/scripts/`](https://github.com/spinalcordtoolbox/PAM50/tree/master/scripts): Copies of up-to-date processing scripts used to generate the PAM50 template and atlas files.
+- [`spinalcordtoolbox/dev/`](https://github.com/spinalcordtoolbox/spinalcordtoolbox-dev/tree/master/dev) ([2016](https://github.com/spinalcordtoolbox/spinalcordtoolbox-dev/tree/265541ec2b1daf121e7a8eb64efd5b764902790c/dev), [2018](https://github.com/spinalcordtoolbox/spinalcordtoolbox-dev/tree/035336b959802f0bb23a3d0b714e674ce5a01e6a/dev)): Archive of older processing scripts used in the general development of SCT (including the PAM50 template and atlas files).

scripts/create_atlas/README.md

@@ -0,0 +1,6 @@
+# To create the atlas of White Matter
+# Run from matlab
+create_atlas.m
+
+# TODO
+- regularize warping fields along z.

scripts/create_atlas/create_masks_csf_and_gm.py

@@ -0,0 +1,85 @@
+#!/usr/bin/env python
+# create masks of CSF and gray matter
+# Author: [email protected]
+# Created: 2014-12-06
+
+# TODO: get GM
+# TODO: add tract corresponding to the undefined values in WM atlas
+
+import glob
+import sys
+import commands
+import numpy as np
+import nibabel as nib
+status, path_sct = commands.getstatusoutput('echo $SCT_DIR')
+# append path that contains scripts, to be able to load modules
+sys.path.append(path_sct + '/scripts')
+import sct_utils as sct
+
+# parameters
+tracts_to_sum_index = 0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29
+folder_atlas = "WMtracts_outputs/final_results"
+file_csf = "WMtract__csf.nii.gz"
+file_gm = "WMtract__gm.nii.gz"
+file_label = 'info_label.txt'
+
+
+def main():
+    # Extract the tracts from the atlas' folder
+    tracts = get_tracts(folder_atlas)
+    nb_tracts = len(tracts)
+    # Get the sum of the tracts
+    tracts_sum = add_tracts(tracts, tracts_to_sum_index)
+    # Save sum of the tracts to niftii
+    save_3D_nparray_nifti(tracts_sum, 'tmp.WM_all.nii.gz', folder_atlas+'/WMtract__00.nii.gz')
+    # binarize it
+    sct.run('fslmaths tmp.WM_all.nii.gz  -thr 0.5 -bin tmp.WM_all_bin.nii.gz')
+    # dilate it
+    sct.run('fslmaths tmp.WM_all_bin.nii.gz  -kernel boxv 5x5x1 -dilM tmp.WM_all_bin_dil.nii.gz')
+    # subtract WM mask to obtain CSF mask
+    sct.run('fslmaths tmp.WM_all_bin_dil -sub tmp.WM_all '+folder_atlas+'/'+file_csf)
+    # add line in info_label.txt
+    text_label = '\n'+str(nb_tracts)+', CSF, '+file_csf
+    open(folder_atlas+'/'+file_label, 'a+b').write(text_label)
+
+def get_tracts(tracts_folder):
+    """Loads tracts in an atlas folder and converts them from .nii.gz format to numpy ndarray 
+    Save path of each tracts
+    Only the tract must be in tracts_format in the folder"""
+    fname_tract = glob.glob(tracts_folder + '/*' + '.nii.gz')
+
+    #Initialise tracts variable as object because there are 4 dimensions
+    tracts = np.empty([len(fname_tract), 1], dtype=object)
+
+    #Load each partial volumes of each tracts
+    for label in range(0, len(fname_tract)):
+       tracts[label, 0] = nib.load(fname_tract[label]).get_data()
+
+    #Reshape tracts if it is the 2D image instead of 3D
+    for label in range(0, len(fname_tract)):
+       if (tracts[label,0]).ndim == 2:
+           tracts[label,0] = tracts[label,0].reshape(int(np.size(tracts[label,0],0)), int(np.size(tracts[label,0],1)),1)
+    return tracts
+
+
+def save_3D_nparray_nifti(np_matrix_3d, output_image, fname_atlas):
+    # Save 3d numpy matrix to niftii image
+    # np_matrix_3d is a 3D numpy ndarray
+    # output_image is the name of the niftii image created, ex: '3D_matrix.nii.gz'
+    img = nib.Nifti1Image(np_matrix_3d, np.eye(4))
+    affine = img.get_affine()
+    np_matrix_3d_nii = nib.Nifti1Image(np_matrix_3d,affine)
+    nib.save(np_matrix_3d_nii, output_image)
+    # copy geometric information
+    sct.run('fslcpgeom '+fname_atlas+' '+output_image, verbose=0)
+
+
+def add_tracts(tracts, tracts_to_sum_index):
+    tracts_sum = np.empty((tracts[0, 0]).shape)
+    for i in tracts_to_sum_index:
+        tracts_sum = tracts_sum + tracts[i, 0]
+    return tracts_sum
+
+
+if __name__ == "__main__":
+    main()

scripts/create_atlas/dnsamplelin.m

@@ -0,0 +1,51 @@
+function out = dnsamplelin(I,n)
+% out = dnsamplelin(I,n)
+% Downsample image I by factor n by computing mean value for each region
+
+I = double(I);
+I = m_normalize(I);
+
+[nx,ny] = size(I);
+nxd = ceil(nx/n);
+nyd = ceil(ny/n);
+out = zeros(nx,ny);
+
+downsamplinggrid = false(size(I));
+downsamplinggrid(1:n:end,1:n:end) = true;
+ind = find(downsamplinggrid);
+[indx,indy] = ind2sub(size(I),ind);
+
+indx = reshape(indx,nxd,nyd);
+indy = reshape(indy,nxd,nyd);
+indx = indx(:,1); indx = indx(:);
+indy = indy(1,:); indy = indy(:);
+xf = indx(end);
+yf = indy(end);
+indx = indx(1:end-1)';
+indy = indy(1:end-1)';
+
+
+for i = indx
+    for j = indy
+        out(i,j) = mean( mean( I(i:i+n-1,j:j+n-1) ) );
+    end
+end
+
+for i = indx
+    out(i,yf) = mean( mean( I(i:i+n-1,yf:end) ) );
+end
+
+for j = indy
+    out(xf,j) = mean( mean( I(xf:end,j:j+n-1) ) );
+end
+
+out(xf,yf) = mean( mean( I(xf:end,yf:end) ) );
+
+out = downsample(out,n);
+out = out';
+out = downsample(out,n);
+out = out';
+
+
+
+end

scripts/create_atlas/m_linear_interp.m

@@ -0,0 +1,21 @@
+function out = m_linear_interp(X,p,q)
+% out = m_linear_interp(X,p,q)
+% Linear interpolation for a 3D data along z (3rd) direction
+% The input data is interpolated between z = n and z = m
+
+X = double(X);
+out = X;
+n = min(p,q);
+m = max(p,q);
+delta = m - n;
+
+if (delta <= 1)
+    out = X;
+else
+    for k = 1:delta-1
+        out(:,:,n+k) = (delta-k)/delta * X(:,:,n) + k/delta * X(:,:,m);
+    end
+end
+
+
+end

scripts/create_atlas/m_normalize.m

@@ -0,0 +1,16 @@
+function out = m_normalize(X)
+% Normalizes input array values [0,1]
+% The class of the output is double
+% If all values in input are equal, they are set to zero
+
+X = double(X);
+
+if( max(X(:)) == min(X(:)) )
+    out = zeros(size(X));
+else
+    out = ( X - min(X(:)) ) / ( max(X(:)) - min(X(:)) );
+end
+
+
+end
+

scripts/create_atlas/m_numbering.m

@@ -0,0 +1,73 @@
+% =========================================================================
+% FUNCTION
+% j_numbering.m
+%
+% Create a cell containing ascending numbers (e.g. '001', '002', ...)
+%
+% INPUTS
+% max_numbering     integer.
+% (nb_digits)       integer. By default, it ajusts it automatically
+% (starting_value)  integer. starting value (default=1)
+% (output_format)   string. 'cell'*, 'array'
+%
+% OUTPUTS
+% out_numbering     cell or array
+%
+% DEPENDANCES
+% (-)
+%
+% COMMENTS
+% Julien Cohen-Adad 2006-11-26
+% =========================================================================
+function varargout = m_numbering(varargin)
+
+
+% initialization
+if (nargin<1) help m_numbering; return; end
+max_numbering = varargin{1};
+if (nargin<2)
+    nb_digits = length(num2str(max_numbering));
+else
+    nb_digits = varargin{2};
+    % check number of digits
+    if (nb_digits<length(num2str(max_numbering)))
+        error('Number of digits too small!!!');
+        return
+    end
+end
+if (nargin<3)
+    starting_value = 1;
+else
+    starting_value = varargin{3};
+end
+if (nargin<4)
+    output_format = 'cell';
+else
+    output_format = varargin{4};
+end
+
+% generate numbering
+out_numbering = cell(max_numbering,1);
+number = starting_value;
+for iNumber=1:max_numbering
+    % write number
+    number_string = num2str(number);
+    % fill with zeros
+    for iDigits=1:nb_digits-length(number_string)
+        number_string = strcat('0',number_string);
+    end
+    out_numbering{iNumber} = number_string;
+    number = number + 1;
+end
+
+if strcmp(output_format,'array')
+    out_numbering_tmp = out_numbering;
+    clear out_numbering
+    for i=1:size(out_numbering_tmp,1)
+		out_numbering(i,:) = out_numbering_tmp{i};
+    end
+	clear out_numbering_tmp
+end
+
+% output
+varargout{1} = out_numbering;

scripts/create_atlas/register_AMU_to_PAM.py

@@ -0,0 +1,113 @@
+#!/usr/bin/env python
+#########################################################################################
+#
+# Register AMU_wm atlas to the PAM_cord segmentation, and also add missing slices at the top and bottom.
+#
+# ---------------------------------------------------------------------------------------
+# Copyright (c) 2016 Polytechnique Montreal <www.neuro.polymtl.ca>
+# Author: [email protected]
+# Created: 2016-08-26
+#
+# About the license: see the file LICENSE.TXT
+#########################################################################################
+
+#TODO: regularize transformations across z
+
+# Import common Python libraries
+import os
+import sys
+import numpy as np
+# append path that contains scripts, to be able to load modules
+path_script = os.path.dirname(__file__)
+path_sct = os.path.abspath(path_script+'/../../../')
+sys.path.append(path_sct+'/scripts')
+import sct_utils as sct
+from msct_image import Image
+
+
+# parameters
+fname_wm = '/Users/julien/Dropbox/Public/sct/PAM50/template/PAM50_wm.nii.gz'
+fname_gm = '/Users/julien/Dropbox/Public/sct/PAM50/template/PAM50_gm.nii.gz'
+fname_cord = '/Users/julien/data/sct_dev/PAM50/template/PAM50_cord.nii.gz'
+
+# create temporary folder
+path_tmp = sct.tmp_create()
+
+# go to temp folder
+os.chdir(path_tmp)
+
+# open volumes
+im_wm = Image(fname_wm)
+data_wm = im_wm.data
+im_gm = Image(fname_gm)
+data_gm = im_gm.data
+im_cord = Image(fname_cord)
+data_cord = im_cord.data
+dim = im_cord.dim
+
+# sum wm/gm
+data_wmgm = data_wm + data_gm
+
+# get min/max z slices from wm/gm
+zsum = np.sum(np.sum(data_wmgm, 0), 0)
+zmin_wm = np.min(np.nonzero(zsum))
+zmax_wm = np.max(np.nonzero(zsum))
+
+# get min/max z slices from cord
+zsum = np.sum(np.sum(data_cord, 0), 0)
+zmin_cord = np.min(np.nonzero(zsum))
+zmax_cord = np.max(np.nonzero(zsum))
+
+# duplicate WM and GM atlas towards the top and bottom slices to match the cord template
+# bottom slices
+for iz in range(zmin_cord, zmin_wm):
+    data_wm[:, :, iz] = data_wm[:, :, zmin_wm]
+    data_gm[:, :, iz] = data_gm[:, :, zmin_wm]
+# top slices
+for iz in range(zmax_wm, zmax_cord):
+    data_wm[:, :, iz] = data_wm[:, :, zmax_wm]
+    data_gm[:, :, iz] = data_gm[:, :, zmax_wm]
+
+# save modified atlases
+im_wm.setFileName('wm_ext.nii.gz')
+im_wm.data = data_wm
+im_wm.save()
+im_gm.setFileName('gm_ext.nii.gz')
+im_gm.data = data_gm
+im_gm.save()
+
+# sum modified wm/gm
+data_wmgm = data_wm + data_gm
+
+# save wm/gm
+im_wm.setFileName('wmgm_ext.nii.gz')
+im_wm.data = data_wmgm
+im_wm.save()
+
+# register wmgm --> cord
+sct.run('cp '+fname_cord+' cord.nii.gz')
+# sct.run('sct_maths -i '+fname_cord+' -laplacian 1 -o cord.nii.gz')
+sct.run('sct_maths -i wmgm_ext.nii.gz -bin 0.5 -o wmgm_ext.nii.gz')
+
+# crop for faster registration
+sct.run('sct_crop_image -i cord.nii.gz -start 40,40,40 -end 100,100,990 -dim 0,1,2 -o cord_crop.nii.gz')
+sct.run('sct_crop_image -i wmgm_ext.nii.gz -start 40,40,40 -end 100,100,990 -dim 0,1,2 -o wmgm_ext_crop.nii.gz')
+# sct.run('sct_maths -i wmgm_ext.nii.gz -laplacian 1 -o wmgm_ext.nii.gz')
+#sct.run('sct_register_multimodal -i wmgm_ext.nii.gz -d cord.nii.gz -iseg wmgm_ext.nii.gz -dseg cord.nii.gz -param step=1,type=im,algo=bsplinesyn,iter=10,slicewise=1,metric=MeanSquares -x linear -r 0')
+sct.run('sct_register_multimodal -i wmgm_ext_crop.nii.gz -d cord_crop.nii.gz -param step=1,type=im,algo=affine,iter=100,slicewise=1,metric=MeanSquares,smooth=1:step=2,type=im,algo=bsplinesyn,iter=5,slicewise=0,metric=MeanSquares,smooth=0 -x linear -r 0')
+sct.run('sct_apply_transfo -i wm_ext.nii.gz -d cord.nii.gz -w warp_wmgm_ext2cord.nii.gz -x linear')
+
+# regularize along S-I direction
+
+# symmetrize
+
+
+# crop below a certain point
+sct.run('sct_crop_image -i wm_ext_reg.nii.gz -dim 2 -start 0 -end 990 -b 0 -o wm_ext_reg_crop.nii.gz')
+
+# rename new file
+sct.run('mv wm_ext_reg_crop.nii.gz PAM50_wm.nii.gz')
+
+# go back to previous folder
+#os.chdir('../')
+

scripts/create_atlas/shift_1pix_to_right.mat

@@ -0,0 +1,4 @@
+1 0 0 -1
+0 1 0 0
+0 0 1 0
+0 0 0 1