Registration tips

[Evelyne]

Hello all!

I was hoping you could give me some tips regarding registration.
My goal is to register fMRI data to the PAM template.
I have acquired a T1 and T2* anatomical image of the spinalcord.

Until now, I have performed these steps:

• Prep T1 image (segment t1, label t1, calculate t1-temp warping field with sct_register_to_template)
• Prep T2 image (segment T2, label t2, calculate t2-temp warping field with sct_register_to_template, calculate GM segmentation)
• Prep fmri images (motion correction and slice time correction).
• Calculate transformation from T2 to T1 using sct multimodal register

Now, I want to register my fmri images to the PAM template.
I am a bit confused about how to do this best.
Initially, I wanted to do something similar to Vahdat et al., 2020: https://doi.org/10.1371/journal.pbio.300078

Then, the warping field needed for the co-registration of the T2*-weighted and T1-weighted images was computed using the SCT registration tools (Fig 1H). Then, by multiplying the warping fields from the T2*-weighted to T1-weighted space (Fig 1H), and from T1-weighted to the template (Fig 1I), we obtained the transformation between the T2*-weighted image and the template. This transformation was further improved by the segmentation of the grey matter in the spinal cord using a dictionary approach [46] and multiplying a corrective local warping field (Fig 1G). This final transformation maps the image from the T2*-weighted space to the template by considering the grey matter structure in the spinal cord. Finally, we computed the parameters needed to co-register the EPI image (Fig 1A) and T2*-weighted image. The resulting transformation (Fig 1F) was then combined with the grey matter corrected transformation from the T2*-weighted space to the template as described earlier, to obtain the final transformation from the EPI to the template (Fig 1E). Note that, in every step, we also calculated an inverse transformation from the destination to the source image, so we also obtained the transformation from the template to the EPI image by multiplying the inverse transformations described earlier.

Is this approach something you would advise?

Looking forward to your tips.

Best,
Evelyne

[joshuacwnewton]

Hi @Evelyne,

Thank you for your question! And, my apologies for the late reply.

SCT has a set of web tutorials, one of which covers the case you describe: Processing fMRI data

Specifically, the page Registering fMRI to the PAM50 template covers how to re-use the existing warping fields from a separate contrast to register fMRI data to the PAM50 template.

((Tagging @jcohenadad to ask for his thoughts about the cited procedure, as well.))

Kind regards,
Joshua

[jcohenadad]

Hi @Evelyne

Thank you for reaching out.

I would advise the approach you are describing only if you are interested in utilizing the GM information from your T2*w image. If not, you can skip the T2*w preprocessing (I assume you meant T2*w preprocessing, and not T2 preprocessing). If you wish to use GM-informed registration (which only works well if the data are of sufficient quality), the last reminding step from what you describe consists in registering the fMRI data to the T2*w data, and then concatenate the transformations fMRI → T2*w → T1w → PAM50.

@joshuacwnewton’s post will help in getting you started with the analysis. Please feel free to reach out with specific questions if you encounter issues/difficulties along the way,

Cheers,
Julien

[Evelyne]

Thank you for your answers.
The first problem I encountered was an error in registering T1 with the template.

This is my code:

if register_t1_temp == 1:
    # define input image (T1)
    input_directory = os.path.join(
        data_pre, subdirectory, "anat")
    wildcard_pattern = 's3*01.nii'
    input_file = fnmatch.filter(os.listdir(
        input_directory), wildcard_pattern)
    input_path = os.path.join(input_directory, input_file[0])

    # define input segmentation
    input_directory_seg = os.path.join(
        data_pre, subdirectory, 'spinalcord', "segment_anat_second")
    wildcard_pattern = '*segmentation.nii'
    input_file_seg = fnmatch.filter(os.listdir(
        input_directory_seg), wildcard_pattern)
    input_path_seg = os.path.join(input_directory_seg, input_file_seg[0])

    # define labels input image
    input_directory_lab = os.path.join(
        data_pre, subdirectory, "spinalcord", "label_anat")
    wildcard_pattern = 's*manlabel.nii'  # name if not manual labelled?
    input_file_lab = fnmatch.filter(os.listdir(
        input_directory_lab), wildcard_pattern)
    input_path_lab = os.path.join(input_directory_lab, input_file_lab[0])

    # define output directory
    output_directory = os.path.join(
        data_pre, subdirectory, "spinalcord", "register_t1_temp")
    # define qc
    qc_directory = os.path.join(output_directory, 'qc_t1_temp')

    # run registering T1 --> template
    command = f"sct_register_to_template -i {input_path} -s {input_path_seg} -l {input_path_lab} -c t1 -qc {qc_directory} -ofolder {output_directory}"
    process = subprocess.Popen(command, shell=True)
    process.wait()
    print('register calc t1 - temp completed for ', subdirectory)

I tested registration with this code:

#test T1 to template
#---------------------------------------------------------------------
if test == 1:
    input_path = '/project/3023009.09/evefra/sub-004/anat/s302300909_sub004-0018-00001-000192-01.nii'
    output_path = '/home/affneu/evefra/sct_5.8/data/PAM50/template/PAM50_t1.nii.gz'
    input_path_warp = '/project/3023009.09/evefra/sub-004/spinalcord/register_t1_temp/warp_anat2template.nii.gz'
    os.chdir('/home/affneu/evefra/')
    command = f"sct_apply_transfo -i {input_path} -w {input_path_warp} -d {output_path}"
    process = subprocess.Popen(command, shell=True)
    process.wait()
    

This is my output:

[jcohenadad]

Hi,

A lot of things could cause this “black image” (eg: issues with the viewer, wrong file selected, issue with path to input data, etc.).

To be able to help, you could send us the necessary data and a SHELL script with relative paths so we can reproduce your analysis and then figure out what went wrong.

[Evelyne]

Hi,

I am pretty sure there is no problem with the viewer, as I tried to visualize it with multiple viewers.
Here is my script:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wednesday 24 may 2023

@author: evefra
"""


# import libraries
import os
import shutil
import subprocess
import fnmatch
import nibabel as nib
import numpy as np
# work in conda env, where SCT is installed. For me: conda activate venv_sct

######################################################################################
# INDICATE WHAT TO DO#
######################################################################################


# ANATOMICAL T1
first_segment_t1 = 0
smooth = 0
second_segment_t1 = 0  # Correct
label_t1 = 0
register_t1_temp = 1  # NOT CORRECT




# LOOP OVER SUB NUMBERS
# ------------------------------------------------------------------------------

sub_list = [20]
for sub in sub_list:

    # data path pre subject unspecific
    data_pre = "/project/3023009.09/evefra"
    subdirectory = f"sub-{sub:03d}"



#################################################################################
        # T1 prep  #
#################################################################################


# FIRST SEGMENTAION OF T1 (including manual C2 and T1 labelling)
# ------------------------------------------------------------------------

    if first_segment_t1 == 1:

        # define input
        input_directory = os.path.join(
            data_pre, subdirectory, "anat")
        wildcard_pattern = 's*.nii'
        input_file = fnmatch.filter(os.listdir(
            input_directory), wildcard_pattern)
        input_path = os.path.join(input_directory, input_file[0])

        # define output
        output_directory = os.path.join(
            data_pre, subdirectory, "spinalcord", "segment_anat")
        if not os.path.exists(output_directory):
            os.makedirs(output_directory)
        output_file = input_file[0].replace("nii", "segmentation.nii")
        output_path = os.path.join(output_directory, output_file)

        # define QC directory
        qc_directory = os.path.join(output_directory, 'qc')

        # run first segmentation
        # label centerline C2 and C7
        command = f"sct_deepseg_sc -i {input_path} -c t1 -o {output_path} -qc {qc_directory} -centerline viewer"
        process = subprocess.Popen(command, shell=True)
        process.wait()
        print('first segment t1 completed for ', subdirectory)


# SMOOTH ALONG SC TO IMPROVE SEGMENTATION
# -------------------------------------------------------------------------------
    if smooth == 1:
        # define input
        input_directory = os.path.join(
            data_pre, subdirectory, "anat")
        wildcard_pattern = 's*.nii'
        input_file = fnmatch.filter(os.listdir(
            input_directory), wildcard_pattern)
        input_path = os.path.join(input_directory, input_file[0])

        # select segmentation image
        input_directory_seg = os.path.join(
            data_pre, subdirectory, 'spinalcord', "segment_anat")
        wildcard_pattern = '*.segmentation.nii'
        input_file_seg = fnmatch.filter(os.listdir(
            input_directory_seg), wildcard_pattern)
        input_path_seg = os.path.join(input_directory_seg, input_file_seg[0])

        # define output
        output_directory = os.path.join(
            data_pre, subdirectory, "spinalcord", "segment_anat")
        output_file_seg = input_file_seg[0].replace(
            "segmentation.nii", "smooth.nii")
        output_path_seg = os.path.join(output_directory, output_file_seg)

        # RUN SMOOOTHING
        command = f"sct_smooth_spinalcord -i {input_path} -o {output_path_seg} -s     {input_path_seg} -smooth 0,0,5"
        process = subprocess.Popen(command, shell=True)
        process.wait()
        print('smooth t1 completed for ', subdirectory)


# DO SECOND SEGMENTATION
# -----------------------------------------------------------------------------

    if second_segment_t1 == 1:

        # define input
        input_directory = os.path.join(
            data_pre, subdirectory, "spinalcord", "segment_anat")
        wildcard_pattern = '*smooth.nii'
        input_file = fnmatch.filter(os.listdir(
            input_directory), wildcard_pattern)
        input_path = os.path.join(input_directory, input_file[0])

        # define output
        output_directory = os.path.join(
            data_pre, subdirectory, "spinalcord", "segment_anat_second")
        if not os.path.exists(output_directory):
            os.makedirs(output_directory)
        output_file = input_file[0].replace(
            "smooth.nii", "second_segmentation.nii")
        output_path = os.path.join(output_directory, output_file)
        # define  qcdirectory
        qc_directory = os.path.join(output_directory, 'qc')

        # select segmentation image
        input_directory_seg = os.path.join(
            data_pre, subdirectory, 'spinalcord', "segment_anat")
        wildcard_pattern = '*.segmentation.nii'
        input_file_seg = fnmatch.filter(os.listdir(
            input_directory_seg), wildcard_pattern)
        input_path_seg = os.path.join(input_directory_seg, input_file_seg[0])

        # run second segmentation
        command = f"sct_deepseg_sc -i {input_path} -c t1 -o {output_path} -qc {qc_directory} -centerline file -file_centerline {input_path_seg}"
        process = subprocess.Popen(command, shell=True)
        process.wait()
        print('second segment t1 completed for ', subdirectory)


# CHECK SEGMENTATION AND CORRECT
# --------------------------------------------------------------------------------

# in FSLeyes




# LABEL T1
# -------------------------------------------------------------------------------
    if label_t1 == 1:

        # define input anatomical image
        input_directory = os.path.join(
            data_pre, subdirectory, "anat")
        wildcard_pattern = 's3*-01.nii'
        input_file = fnmatch.filter(os.listdir(
            input_directory), wildcard_pattern)
        input_path = os.path.join(input_directory, input_file[0])

        # define input segmentation
        input_directory_seg = os.path.join(
            data_pre, subdirectory, 'spinalcord', "segment_anat_second")
        wildcard_pattern = '*segmentation.nii'
        input_file_seg = fnmatch.filter(os.listdir(
            input_directory_seg), wildcard_pattern)
        input_path_seg = os.path.join(input_directory_seg, input_file_seg[0])

        # define output paths and name
        output_directory = os.path.join(
            data_pre, subdirectory, "spinalcord", "label_anat")
        if not os.path.exists(output_directory):
            os.makedirs(output_directory)
        output_file_manual = input_file[0].replace(".nii", "manlabel.nii")
        output_path_manual = os.path.join(output_directory, output_file_manual)
        output_file = input_file[0].replace(".nii", "label.nii")
        output_path = os.path.join(output_directory, output_file)

        # define qc directory
        qc_directory_man = os.path.join(output_directory, 'qc_manual')
        qc_directory = os.path.join(output_directory, 'qc')

        # manual label C2 and C7
        command = f"sct_label_utils -i {input_path} -qc {qc_directory_man} -create-viewer 3,7 -o {output_path_manual} -msg 'Please manually label the inter-vertebral disc '"
        process = subprocess.Popen(command, shell=True)
        process.wait()

    # define input man label
    # wildcard_pattern = '*manlabel.nii'
    # input_file_man = fnmatch.filter(os.listdir(
    #   output_directory), wildcard_pattern)
    # input_path_man = os.path.join(output_directory, input_file_man[0])

    # correcting s or qform if needed
    # sct_image -set-sform-to-qform -i {

    # run labelling
    # command = f"sct_label_vertebrae -i {input_path} -c t1 -s {input_path_seg} -qc {qc_directory} -o {output_path} -initlabel {input_path_man}"
    # process = subprocess.Popen(command, shell=True)
    # process.wait()
        print('label t1 completed for ', subdirectory)



# REGISTER FROM T1 to TEMPLATE
# -----------------------------------------------------------------------------

    if register_t1_temp == 1:
        # define input image (T1)
        input_directory = os.path.join(
            data_pre, subdirectory, "anat")
        wildcard_pattern = 's3*01.nii'
        input_file = fnmatch.filter(os.listdir(
            input_directory), wildcard_pattern)
        input_path = os.path.join(input_directory, input_file[0])

        # define input segmentation
        input_directory_seg = os.path.join(
            data_pre, subdirectory, 'spinalcord', "segment_anat_second")
        wildcard_pattern = '*segmentation.nii'
        input_file_seg = fnmatch.filter(os.listdir(
            input_directory_seg), wildcard_pattern)
        input_path_seg = os.path.join(input_directory_seg, input_file_seg[0])

        # define labels input image
        input_directory_lab = os.path.join(
            data_pre, subdirectory, "spinalcord", "label_anat")
        wildcard_pattern = 's*manlabel.nii'  # name if not manual labelled?
        input_file_lab = fnmatch.filter(os.listdir(
            input_directory_lab), wildcard_pattern)
        input_path_lab = os.path.join(input_directory_lab, input_file_lab[0])

        # define output directory
        output_directory = os.path.join(
            data_pre, subdirectory, "spinalcord", "register_t1_temp")
        # define qc
        qc_directory = os.path.join(output_directory, 'qc_t1_temp')

        # run registering T1 --> template
        command = f"sct_register_to_template -i {input_path} -s {input_path_seg} -l {input_path_lab} -c t1 -qc {qc_directory} -ofolder {output_directory}"
        process = subprocess.Popen(command, shell=True)-
        process.wait()
        print('register calc t1 - temp completed for ', subdirectory)

I will send the data private to you.
Thank you in advance.
Best,
Evelyne

[jcohenadad]

Hi @Evelyne,

I tried putting your script in a file and running it (after declaring SCT’s conda environment), however I ran into an issue:

  File "/Users/julien/Desktop/share/script.py", line 255
    process = subprocess.Popen(command, shell=True)-
                                                    ^
SyntaxError: invalid syntax

Moreover, I noticed that your script includes hard-coded paths that don’t match my environment, eg:

    data_pre = "/project/3023009.09/evefra"

I do not have the time to dig in your script and adapt it to my local environment, so could you please send me a script that works “out of the box” with the files you sent me in private?

On another note, I opened the files you sent me and noticed that the image and the segmentation are not in the same physical space, see the different qform from the NIfTI headers:

fslhd s302300909_sub020-0018-00001-000192-01.nii
...
qto_xyz:1	0.000000 0.000000 1.000000 -95.500000 
qto_xyz:2	-1.000000 0.000000 -0.000000 146.947952 
qto_xyz:3	0.000000 1.000000 -0.000000 -63.840141 
qto_xyz:4	0.000000 0.000000 0.000000 1.000000 

versus:

fslhd s302300909_sub020-0018-00001-000192-01.segmentation.nii
...
qto_xyz:1	0.000000 0.000000 1.000000 -94.804344 
qto_xyz:2	-1.000000 0.000000 -0.000000 116.308746 
qto_xyz:3	0.000000 1.000000 -0.000000 -174.108688 
qto_xyz:4	0.000000 0.000000 0.000000 1.000000 

After digging a bit more, I noticed a mismatch between the sform and the qform of your anatomical image:

qform_name	Scanner Anat
qform_code	1
qto_xyz:1	0.000000 0.000000 1.000000 -95.500000 
qto_xyz:2	-1.000000 0.000000 -0.000000 146.947952 
qto_xyz:3	0.000000 1.000000 -0.000000 -63.840141 
qto_xyz:4	0.000000 0.000000 0.000000 1.000000 
qform_xorient	Anterior-to-Posterior
qform_yorient	Inferior-to-Superior
qform_zorient	Left-to-Right
sform_name	Scanner Anat
sform_code	2
sto_xyz:1	0.000000 0.000000 1.000000 -94.804344 
sto_xyz:2	-1.000000 0.000000 0.000000 116.308746 
sto_xyz:3	0.000000 1.000000 0.000000 -174.108688 
sto_xyz:4	0.000000 0.000000 0.000000 1.000000 
sform_xorient	Anterior-to-Posterior
sform_yorient	Inferior-to-Superior

This is the likely cause of your issues. In order to solve this, you could copy the sform into the qform and then re-run the segmentation:

# Set sform to qform
sct_image -i s302300909_sub020-0018-00001-000192-01.nii -set-sform-to-qform -o s302300909_sub020-0018-00001-000192-01_qform.nii
# Run segmentation
sct_deepseg_sc -i s302300909_sub020-0018-00001-000192-01_qform.nii -c t1 -qc qc

In this case the output qform of the image and of the segmentation do match:

fslhd s302300909_sub020-0018-00001-000192-01.nii 
...
qto_xyz:1	0.000000 0.000000 1.000000 -95.500000 
qto_xyz:2	-1.000000 0.000000 -0.000000 146.947952 
qto_xyz:3	0.000000 1.000000 -0.000000 -63.840141 
qto_xyz:4	0.000000 0.000000 0.000000 1.000000 

fslhd s302300909_sub020-0018-00001-000192-01_qform_seg.nii
...
qto_xyz:1	0.000000 0.000000 1.000000 -95.500000 
qto_xyz:2	-1.000000 0.000000 -0.000000 146.947952 
qto_xyz:3	0.000000 1.000000 -0.000000 -63.840141 
qto_xyz:4	0.000000 0.000000 0.000000 1.000000 

Let me know if that solves your issue.

[joshuacwnewton]

Hi @Evelyne,

I am another SCT developer, and I’m helping @jcohenadad to debug this issue.

@jcohenadad

I do not have the time to dig in your script and adapt it to my local environment, so could you please send me a script that works “out of the box” with the files you sent me in private?

I tried out @Evelyne’s script, and the necessary fixes are quite simple:

1. Rearrange the files to have the following structure:
   
2. Change data_pre to point to the parent folder for the newly-created sub-020 folder.
3. Remove the - character that is causing the syntax error.

After doing this, I am indeed able to run the script and reproduce the black image issue.

---

This is the likely cause of your issues. In order to solve this, you could copy the sform into the qform and then re-run the segmentation.

I tried your suggestion (though I had to format the output name differently to match the expectations of the script):

sct_image -i s302300909_sub020-0018-00001-000192-01.nii -set-sform-to-qform \
          -o s302300909_qform_sub020-0018-00001-000192-01.nii

Then, I deleted the previously-generated intermediate files (everything but anat), and re-ran the script with various settings turned on (to regenerate the segmentation):

# ANATOMICAL T1
first_segment_t1 = 1
smooth = 1
second_segment_t1 = 1  # Correct
label_t1 = 1
register_t1_temp = 1  # NOT CORRECT

This requires some manual input for sct_deepseg_sc -centerline-viewer and the disc 3/7 labeling, so I tried to recreate the labels as faithfully as possible.

With the qform/sform fix, I can confirm that the PAM50 template is warped to the subject image and vice-versa. However, the quality of the registration is somewhat poor (mismatched vertebrae), so the registration parameters will likely need to be tweaked:

I will now try to investigate why this qform/sform mismatch causes the registration failure in the first place.

Kind regards,
Joshua

[Evelyne]

Thankyou. It works now.

[Evelyne]

When I try to concatenate the warping fields from t1-temp and t2-t1, I come across the following problem: The dimensions are not the same. I already tried to cut the images with no success. I have now overcome this issue by not doing the concatenation, but just selecting the three warping fields in the sct_apply_transfo function. However, this restricts me from adding grey matter segmentation to improve registration. Any thoughts on this? This is my code and two examples:

    concatenated_warping_field = '/project/3023009.09/evefra/sub-024/spinalcord//project/3023009.09/evefra/sub-024/spinalcord/concat_t2_t1_temp/warp.nii.gz'
    
command = f"sct_image -i /project/3023009.09/evefra/sub-024/spinalcord/register_t1_temp/warp_anat2template.nii.gz /project/3023009.09/evefra/sub-024/spinalcord/register_t2_t1/warp_s302300909_sub024-0009-00001-000001-012s302300909_sub024-0016-00001-000192-01.qform.nii.gz'
    -concat t -o {concatenated_warping_field}"
    process = subprocess.Popen(command, shell=True)
    process.wait()

This is the error:

--
Spinal Cord Toolbox (5.8)

sct_image -i /project/3023009.09/evefra/sub-024/spinalcord/register_t1_temp/warp_anat2template.nii.gz /home/affneu/evefra/Spinalcord/scripts/cropped_image.nii.gz -concat t -o /project/3023009.09/evefra/sub-024/spinalcord//project/3023009.09/evefra/sub-024/spinalcord/concat_t2_t1_temp/warp.nii.gz
--

Traceback (most recent call last):
  File "/home/affneu/evefra/sct_5.8/spinalcordtoolbox/scripts/sct_image.py", line 610, in <module>
    main(sys.argv[1:])
  File "/home/affneu/evefra/sct_5.8/spinalcordtoolbox/scripts/sct_image.py", line 254, in main
    im_out = [concat_data(im_in_list, dim)]
  File "/home/affneu/evefra/sct_5.8/spinalcordtoolbox/image.py", line 909, in concat_data
    data_concat = np.concatenate(dat_list, axis=dim)
  File "<__array_function__ internals>", line 180, in concatenate
ValueError: all the input array dimensions for the concatenation axis must match exactly, but along dimension 2, the array at index 0 has a size of 991 and the array at index 1 has a size 192

Data: SURFdrive

[Evelyne]

Is it for example possible to select multiple -initwarp in the sct_register_multimodal function?

[joshuacwnewton]

Hi @Evelyne,

My deepest apologies for the late reply.

@Evelyne

Is it for example possible to select multiple -initwarp in the sct_register_multimodal function?

I don’t believe this is currently possible. I think you are taking the correct approach (by concatenating the warping fields).

@Evelyne

This is my code and two examples:

sct_image -i /project/3023009.09/evefra/sub-024/spinalcord/register_t1_temp/warp_anat2template.nii.gz /home/affneu/evefra/Spinalcord/scripts/cropped_image.nii.gz -concat t -o /project/3023009.09/evefra/sub-024/spinalcord//project/3023009.09/evefra/sub-024/spinalcord/concat_t2_t1_temp/warp.nii.gz

Instead of using the sct_image -concat t function, I would recommend using the sct_concat_transfo function instead. ((Warping fields are actually 5D images (as described on this page of SCT’s documentation), so regular image concatenation via the t dimension will not work in this case. Hence, why SCT provides a specialized script for concatenating warping fields.))

Kind regards,
Joshua

[Evelyne]

Thankyou @joshuacwnewton

Did you already figure out why a blank screen occurred? I get another blank screen when improving my registration using sct_multimodal_transfo and gm segmentations. Without the addition of gm, the warping field is working. Some data (t2, initial warp, my result): SURFdrive

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wednesday 24 may 2023

@author: evefra
"""


# import libraries
import os
import shutil
import subprocess
import fnmatch
import nibabel as nib
import numpy as np
# work in conda env, where SCT is installed. For me: conda activate venv_sct

######################################################################################
# INDICATE WHAT TO DO#
######################################################################################



# COREGISTERING
register_t2_t1 = 0
register_fmri_t2 = 0
t2_to_t1_mult_t1_to_temp = 0
correct_qform_sc = 0
gm_segment = 0
gm_segment_add = 0
epi_to_t2_mult_t2_to_temp = 0
final_warp_apply = 0

# LOOP OVER SUB NUMBERS
# ------------------------------------------------------------------------------

sub_list = [10]
# [1, 2, 4, 5, 6, 8, 9, 10, 11, 13, 14, 16, 17, 18, 24]
for sub in sub_list:

# CORRECT QFORM ANATOMICAL
# -----------------------------------------------------------------------------

    if correct_qform_sc == 1:

        # Define the input
        input_directory_d = os.path.join(
            data_pre, subdirectory, "anat_sc")
        wildcard_pattern_d = 's3*01.nii'
        input_file_d = fnmatch.filter(os.listdir(
            input_directory_d), wildcard_pattern_d)
        input_path_des = os.path.join(input_directory_d, input_file_d[0])

        # define output
        output_directory = os.path.join(
            data_pre, subdirectory, "anat_sc")
        output_file = input_file_d[0].replace("nii", "qform.nii")
        output_path = os.path.join(output_directory, output_file)

        # run command
        command = f"sct_image -i {input_path_des} -set-sform-to-qform -o {output_path}"
        process = subprocess.Popen(command, shell=True)
        process.wait()

        print('qform fix anat sc completed for ', subdirectory)
        
# CALC SEGMENTATION GM
# ----------------------------------------------------------------------------

    if gm_segment == 1:

        # define input
        input_directory_d = os.path.join(
            data_pre, subdirectory, "anat_sc")
        wildcard_pattern_d = 's3*01.qform.nii'
        input_file_d = fnmatch.filter(os.listdir(
            input_directory_d), wildcard_pattern_d)
        input_path_des = os.path.join(input_directory_d, input_file_d[0])
       
        # define output
        output_directory = os.path.join(
            data_pre, subdirectory, "spinalcord", "gm_segment_t2_anat")
        output_file_seg = 'gm_t2.nii'
        output_path = os.path.join(output_directory, output_file_seg)

        # Check if the removed volumes directory exists, create it if necessary
        if not os.path.exists(output_directory):
            os.makedirs(output_directory)

        # define qc
        qc_directory = os.path.join(output_directory, 'qc')

        # Run gm segment
        command = f"sct_deepseg_gm -i {input_path_des} -qc {qc_directory} -o {output_path}"
        process = subprocess.Popen(command, shell=True)
        process.wait()
        print('segment gm completed for ', subdirectory)


# ADD SEGMENTATION GM
# ----------------------------------------------------------------------------
    if gm_segment_add == 1:

        # define input template (always the same)
        input_directory = '/home/affneu/evefra/sct_5.8/data/PAM50/template/PAM50_t1.nii.gz'

        # define input gm segmentation (always the same)
        input_directory_iseg = '/home/affneu/evefra/sct_5.8/data/PAM50/template/PAM50_gm.nii.gz'

        # define destination (T2)
        input_directory_d = os.path.join(
            data_pre, subdirectory, "anat_sc")
        wildcard_pattern_d = 's3*01.nii'
        input_file_d = fnmatch.filter(os.listdir(
            input_directory_d), wildcard_pattern_d)
        input_path_des = os.path.join(input_directory_d, input_file_d[0])

        # define destination segmentation (T2 GM)
        input_directory_des_seg = os.path.join(
            data_pre, subdirectory, "spinalcord", "gm_segment_t2_anat")
        wildcard_pattern_des_seg = 'gm*.nii'
        input_file_des_seg = fnmatch.filter(os.listdir(
            input_directory_des_seg), wildcard_pattern_des_seg)
        input_path_des_seg = os.path.join(
            input_directory_des_seg, input_file_des_seg[0])

        # define init warp
        input_directory_warp = os.path.join(
            data_pre, subdirectory, "spinalcord", "concat_t2_t1_temp")
        wildcard_pattern_warp = 'warp_t22temp.nii.gz'
        input_file_warp = fnmatch.filter(os.listdir(
            input_directory_warp), wildcard_pattern_warp)
        input_path_des_warp = os.path.join(
            input_directory_warp, input_file_warp[0])

        # define output warp
        output_directory = os.path.join(
            data_pre, subdirectory, "spinalcord", "concat_t2_t1_temp")
        if not os.path.exists(output_directory):
            os.makedirs(output_directory)
        output_file = 'gm_warp_t22temp.nii.gz'
        output_path_warp = os.path.join(output_directory, output_file)
        
        print(input_directory)
        print(input_directory_iseg)
        print(input_path_des)
        print(input_path_des_seg)
        print(input_path_des_warp)
        print(output_path_warp)        

        # run warping calculation
        command = f"sct_register_multimodal -i {input_directory} -iseg {input_directory_iseg} -d {input_path_des} -dseg {input_path_des_seg} -initwarp {input_path_des_warp} -owarp {output_path_warp}"
        process = subprocess.Popen(command, shell=True)
        process.wait()
        print('add gm to t2-temp warp completed for ', subdirectory)

[joshuacwnewton]

@Evelyne

Did you already figure out why a blank screen occurred?

Yes, we did! The earlier blank image was caused by the mismatched qform and sform matrices in the input image. If you are interested, there is a more technical description of the internal cause of the error here. But, in short, the best way to avoid this is by making sure that the sform and qform matrices match in your input anat image before producing the segmentation/labels/registration/etc.

However, looking at the sample images you shared, it looks like you already solved the sform/qform issues. So, I will begin investigating the new issue with sct_register_multimodal now.

Kind regards,
Joshua

[joshuacwnewton]

Hi @Evelyne,

It looks like the script you shared in this comment is only partially copied and pasted. Some of the steps listed below are missing:

# COREGISTERING
register_t2_t1 = 0
register_fmri_t2 = 0
t2_to_t1_mult_t1_to_temp = 0
correct_qform_sc = 0
gm_segment = 0
gm_segment_add = 0
epi_to_t2_mult_t2_to_temp = 0
final_warp_apply = 0

Could you please share the entire processing script by uploading your “.py” file in your reply? (I’d like to try running the full pipeline (segmentation, labeling, registration, etc.) from start to finish, just to make sure I can reproduce the issue on my machine.)

Kind regards,
Joshua

[Evelyne]

Dear @joshuacwnewton,

Thank you for helping out.
Hereby the script:
sct_preproces.py (37.7 KB)

[joshuacwnewton]

@Evelyne

I get another blank screen when improving my registration using sct_multimodal_transfo and gm segmentations. Without the addition of gm, the warping field is working.

I tried out your script with “gm_segment = 1” and “gm_segment_add = 1” plus the input files you provided. I did have to make 1 small change to your script, though:

    if gm_segment_add == 1:

        # define destination (T2)
        input_directory_d = os.path.join(
            data_pre, subdirectory, "anat_sc")
-        wildcard_pattern_d = 's3*01.nii'
+        wildcard_pattern_d = 's3*01.qform.nii'

After making this change, I got the following results:

• I got an empty/blank “PAM50_T1_reg.nii.gz” file, just as you describe.
• The GM-informed warping field “gm_warp_t22temp.nii.gz” is also empty (all voxels are 3.4028234663852886e+38.)
• The input GM segmentation looks good.
• However, the concatenated -initwarp warping field “warp_t22temp.nii.gz” is empty (all voxels are 3.4028234663852886e+38).

So, I think the “sct_concat_transfo” command (that produced the file “warp_t22temp.nii.gz”) might be the issue here. However, I’m unable to troubleshoot the “t2_to_t1_mult_t1_to_temp” step in your script because I don’t have the necessary files to reproduce your steps. I was mistaken. Please see this comment.

Could you please provide all of the input files for sub-010 (both T1 and T2, and anything else that may be necessary, if possible) so that I can follow each step in the script from start to finish?

Thank you kindly,
Joshua

[Evelyne]

Dear @joshuacwnewton

Thank you for helping out.
I have also uploaded a T1 and one EPI image.

I do not think the problem is in warp_t22temp.nii.gz. Whenever I apply this warp, I can successfully register the T2 to PAM50. I think the problem is in the next step when the initial warp is combined with the GM segmentation.

[joshuacwnewton]

@Evelyne

I do not think the problem is in warp_t22temp.nii.gz. Whenever I apply this warp, I can successfully register the T2 to PAM50. I think the problem is in the next step when the initial warp is combined with the GM segmentation.

My deepest apologies! You are absolutely correct. I had forgotten to adjust my contrast settings on my image viewer (FSLeyes), and so warp_t22temp.nii.gz appeared empty, when it was indeed not.

Looking into this further now. I will reply when I have a further update.

[joshuacwnewton]

Hi @Evelyne,

I took a closer look at the final sct_register_multimodal command inside the “gm_segment_add” step, and compared it to the command used in the GM-informed registration between the PAM50 template and T2* data - Spinal Cord Toolbox documentation page of the SCT tutorials.

I noticed the following differences:

1. For GM-informed registration, the WM segmentation should be used (to retain the shape of both the cord and the GM)
2. The version of the PAM50 template used should match the contrast of the -d image (in this case, I think PAM50_t2.nii.gz should be used instead of PAM50_t1.nii.gz, specifically because we are using a T2 anat image here).
3. A different set of registration parameters are used for -param, due to this being a “fine-tuning” registration thanks to the use of -initwarp (rather than a more coarse registration).
4. In sct_register_multimodal, -initwarp gets applied to the source image (-i), and -initwarpinv gets applied to the dest image (-d).
   • In your command, warp_t22temp.nii.gz (i.e. anat->template) is passed to -initwarp, so it will be applied to the source image.
   • However, the source image (-i) in this case is the template image. So, I believe the wrong initial warping field is being applied to the PAM50 template.

I fixed the 1st, 2nd, and 3rd discrepancies. However, this still resulted in a blank registration warping field.

I then fixed the 4th discrepancy by swapping the source and dest images, that way warp_t22temp.nii.gz gets applied to the T2 image and not the PAM50 template image. This successfully generated a new GM-informed gm_warp_t22temp.nii.gz:

However, one small detail to mention regarding my fix is that, because only -initwarp was provided, only the anat->template warping field will be generated. If you want to also produce the inverse (template->anat) warping field, you will need to supply a second initial warping field (i.e. warp_temp2t2.nii.gz), which I think will mean going back and adding a second sct_concat_transfo command to your earlier step.

Anyway, here are the total changes I made to your script: Changes made to sct_preprocess.py - Diff Checker (Note: Please disregard the filepath and subprocess.Popen changes I made! I needed to make those changes to get the script to work on my machine.)

Kind regards,
Joshua