Lesion segmentation analysis

[plbenveniste]

In this issue, I use the segmentation of the lesions on the spinal cord after cropping the images (issue #37 ).

The analysis of the segmentations accounts for the following:

• Identifying a lesion across slides
• Measuring its volume
• Measuring its center
• Registration between images and segmentations at different time points
• Comparison of the lesions at different time points
• Creating a text file to sum up evolutions of the MS lesions

[plbenveniste]

I worked on a segmentation that was inferred by nnU-Net (the one trained on 1000 epochs) on image sub-cal072_ses-M12_STIR. On this volume, I performed:

• Automatic clustering of the segmented voxels in lesions (the number of lesion is automatically detected) across slices
• Computation of the lesion volumes
• Computation of the center of the lesions

Output:

It outputs a file where each segmented voxel is indexed per lesion (1 for lesion 1, 2 for lesion 2, etc.)

Number of lesions: 4
Volume and center of each lesion (mm3):
Lesion 1 : volume: 42.63 mm3, center: [155.5862069 224.79310345 4. ]
Lesion 2 : volume: 232.26 mm3, center: [159.92405063 170.44303797 4.34810127]
Lesion 3 : volume: 29.4 mm3, center: [164.4 132.3 3. ]
Lesion 4 : volume: 29.4 mm3, center: [218.75 21.5 4. ]

[valosekj]

• Automatic clustering of the segmented voxels in lesions (the number of lesion is automatically detected) across slices
• Computation of the lesion volumes
• Computation of the center of the lesions

Cool! For this, you are using the lesion_seg_analysis.py script from plb/nnunet branch, right?

Be aware that SCT contains sct_analyze_lesion.py function, which does pretty similar things:

$ sct_analyze_lesion -m sub-cal080_ses-M0_STIR_lesion-manual.nii.gz 

...

Lesion count = 1

Measures on lesion #1...
  Volume : 52.92 mm^3

...

Averaged measures...
  Volume [mm^3] = 52.92+/-0.0
  (S-I) Length [mm] = nan+/-nan
  Equivalent Diameter [mm] = nan+/-nan

Total volume = 52.92 mm^3
Lesion count = 1

...

(sct_analyze_lesion.py does not print the lesion center, but this could be easily added)

[plbenveniste]

Currently performing registration between M0 and M12 images using:
sct_register_multimodal -i sub-cal072_ses-M0_STIR.nii.gz -d sub-cal072_ses-M12_STIR.nii.gz -iseg sub-cal072_ses-M0_STIR_seg.nii.gz -dseg sub-cal072_ses-M12_STIR_seg.nii.gz with the segmentations being the spinal cord segmentation (computed with sct_deepseg_sc) used to improve the registration.

The registration is not of great quality because of the very limited number of slices on the z-axis (9 slices in this case). This will surely have a clear impact on the study of the evolution of spinal cord lesions.

[jcohenadad]

what is the purpose of the registration? if it is to identify the same lesions across time points, then it should be fine (given that lesions span several voxels)

[plbenveniste]

The objective is to register images from M0 and M12 to perform comparison of the sc lesions. Yes, it means identifying the lesions across time points.

[valosekj]

Currently performing registration between M0 and M12 images using: sct_register_multimodal -i sub-cal072_ses-M0_STIR.nii.gz -d sub-cal072_ses-M12_STIR.nii.gz -iseg sub-cal072_ses-M0_STIR_seg.nii.gz -dseg sub-cal072_ses-M12_STIR_seg.nii.gz with the segmentations being the spinal cord segmentation (computed with sct_deepseg_sc) used to improve the registration.

Another possible approach might be using sct_straighten_spinalcord; see point 2 in neuropoly/idea-projects#17 (comment). However, sct_straighten_spinalcord requires, besides SC seg, also disc labels (which are not currently available for STIR/PSIR).

[jcohenadad]

from @valosekj:

Another possible approach might be using sct_straighten_spinalcord; see point 2 in neuropoly/idea-projects#17 (comment). However, sct_straighten_spinalcord requires, besides SC seg, also disc labels (which are not currently available for STIR/PSIR).

That's a good idea, although it might be an overkill (straightening takes time). But for visualisation purpose, that could be interesting. In any case, if @plbenveniste needs to have some reasonable level of alignment between time points, he will definitely need disc levels (otherwise, how can he make sure that the coverage is exactly the same across time points).

from @plbenveniste:

The objective is to register images from M0 and M12 to perform comparison of the sc lesions. Yes, it means identifying the lesions across time points.

In that case, you don't need sub-millimetric precision, given that lesions are 'big blobs'. 2-3mm precision in the registration will make it possible to identify the same lesion across time.

[plbenveniste]

After registration, there is still one problem. Even though registration is well performed on slices, because of the anisotropic nature of the image, we can see slices with spinal cord "registered" to slices without spinal cord. That happens, even though both images are well aligned on the x and y axis.

The resulting problem is visible when performing registration on the lesion segmentation images. The time points can't really be compared since they don't superpose well.

Left image = 1st time point image / Right image = 2nd time point image registered to 1st time point
Red labels = 1st time point lesions / Blue labels = 2nd time point labels (registered to 1st time point)

Solution:

• One solution is to perform lesion comparison from M0 to M12 without doing registration, hoping that lesions from the first time point are "similar enough" to those at the second time point and close enough.
• Another solution I want to try is to create a .nii.gz image in which 0=spinal cord and 1=lesion segmentation. From this, I'll directly perform registration to superpose lesions and then compare them.

[jcohenadad]

I'm not sure I fully understand #38 (comment). A few comments:

• you mentioned "the registration is well performed on slices", but the first GIF animation suggests that the two time points are not well registered (notice a shift of ~5cm along the rostro-caudal axis)
• also, what do you mean by registration performed "on slices"? How was the registration done by the way (linking to your code would be useful)
• "both images are well aligned on the x and y axis": in medical imaging, I would recommend to instead use the terms AP, RL, IS to avoid confusion (because an "x" axis can be along AP, RL or IS depending on the orientation of the FOV during acquisition)
• can you make the images brighter/more exposure? it is very difficult to properly see these images
• what are the red and blue labels on the GIF images?

[plbenveniste]

The code is here

[jcohenadad]

additional comments:

• always use QC report
• make sure to add params in your registration, notably: use the seg as step1 (with centermass or slicereg-- play with polynom order), with syn (step2) make sure to use deformation along 1x1x1, metrics: should NOT be meansquare (bc images aquired at diff time points) but instead CC, experience with smoothing (useful in noisy datasets), make sure to use -x linear (or sync)
• you could potentially use the vert labels (but start without-- i think you should be fine without)
• you might be able to register without the segmentation, using the "dl" algo, and if that doesn't work, start with translation and also try syn with various shrink factor.

[plbenveniste]

Changing the parameters improved the registration. The final parameters are:

step=1,type=seg,algo=slicereg,poly=5:step=2,type=im,algo=syn,metric=MI,deformation=1x1x1,smooth=3

I played with the polynomial order, the metric (MI or CC), and the smoothing.

However, even though the spinal cord is well aligned, there is still a problem with the alignment of vertebral labels as we can see in the following:

We can see that the first time point is not overlapping correctly with the second time point, especially when looking at the vertebral levels (look at the green cursor) or when looking at the label from the first time point on both the image from the first time point and the image from the second time point.

I am now trying to add vertebral levels to improve the quality.

[plbenveniste]

I can't manage to use vertebral labels in the following command line.

The error I get is the following:

[valosekj]

Exception: Error: number of source and destination landmarks are not the same, so landmarks cannot be paired.

Do the -ilabel and -dlabel have the same number of labels?

[plbenveniste]

Thanks for the tip. One of the vertebral label files had one more label.

Now I only keep common labels in both files before using them for registration.

The registration is much better when using vertebral levels. We can see the lesions clearly overlapping each other.

Now, moving on to lesion comparison and temporal evolution of the lesions.

[jcohenadad]

now that looks like a much better registration indeed!

[jcohenadad]

can you also try with the 'dl' algo (without segmentation)? I'm curious to see how it performs

[plbenveniste]

Registration performed with 'syn' algo without the spinal cord segmentation but with vertebral levels:

parameters = 'step=0,type=label,dof=Tx_Ty_Tz:step=1,type=im,algo=syn,metric=CC,deformation=1x1x1,smooth=3'

Registration performed with 'dl' with vertebral levels:

parameters = 'step=0,type=label,dof=Tx_Ty_Tz:step=1,type=im,algo=dl,metric=CC,deformation=1x1x1,smooth=3'

Registration performed with 'dl' without the vertebral levels:

parameters = 'step=1,type=im,algo=dl,metric=CC,deformation=1x1x1,smooth=3'

Conclusion :

It seems that overall the best results are obtained with the 'dl' algorithm with vertebral levels. Even though, it doesn't require spinal cord segmentation, computing vertebral levels does: so in terms of computing time, both 'syn' and 'dl' seem similar.

Next objectives:

• Registration is done with algo = 'dl' and uses vert levels
• Identify lesions on both time points
• Create the corresponding matching between lesions through time
• Compute the size and the centers of each lesion for time point 0
• Change the label masks: the value of the voxel = ID of the lesion for the
• Apply the inverse warping field on the lesion mask of the second time point
• Compute the volume and center of each lesion
• Generate a '.txt' file with the output of the analysis

[valosekj]

parameters = 'step=0,type=label,dof=Tx_Ty_Tz:step=1,type=im,algo=dl,metric=CC,deformation=1x1x1,smooth=3'

I believe that when algo=dl is used, all other parameters (such as metric=CC,deformation=1x1x1,smooth=3) are ignored. --> There is no need to specify them. -param step=1,type=im,algo=dl should be enough.

[plbenveniste]

Registration, lesion segmentation, and lesion matching (from first-time point to 2nd-time point) have been done. I did these using the lesions from the 2nd time point on which forward had been applied. Everything is therefore viewed on the image from the first time point. I created a lesion segmentation file where the voxel value = lesion ID. The objective is to be able to use a matching table to match lesions from the 2nd time point to the 1st. I then applied inverse warping on this lesion mask. The result is the following:

The inverse warping also modified the intensity of the voxels which corrupts the labeling of the lesions as IDs.

In order to avoid that problem, I performed lesion clustering on the first time point and the 2nd time point warped on the first time point image. I got the temporal matching of lesions this way.

Then, I did clustering of the lesions in the second time point without warping (meaning that they are still in their original image). To get the matching between the second time point images with and without warping, I used graph-based matching. I created a graph where each node is the center of the lesion and each edge is the distance between the lesions. I, therefore, obtained the matching without having to perform inverse warping.

TO DO:

• include the volume of the lesion in the node embedding.
• generate the final '.txt' output

[jcohenadad]

The inverse warping also modified the intensity of the voxels which corrupts the labeling of the lesions as IDs.

with nearest neighbour interpolation, the intensity of the voxel should remain the same

[plbenveniste]

After quite some trouble with graph matching, I changed my strategy to the following:

• perform registration between ses-1 and ses-2: the objective is to get the segmentation file from ses-2 registered to the segmentation file from ses-1
• perform lesion matching on ses-1 and ses-2_registered_to_ses-1: both segmentation file have voxels = lesion ID (with IDs matching from ses-1 to ses-2)
• perform inverse warping on ses-2_registered_to_ses-1_correctly_labeled
• perform lesion matching between the original ses-2 segmentation and the inversely wrapped segmentation of ses-2: the original ses-2 seg file is transformed so that each voxel = lesion ID (with IDs matching from ses-1 to ses-2)
• perform a comparison of the lesions from ses-1 to ses-2
• create an output file

The result is the following:

On the left the ses-1 image, on the right the ses-2 image (both with matching labeled lesions)

The output file is the following:

The code can be found here

[jcohenadad]

very good progress @plbenveniste ! 👏

On the left the t1 image, on the right the t2 image

please avoid using t1 t2 for times points 1-2, because in MRI we often use t1 and t2 for the contrasts. Instead, you could use ses-1, ses-2 (inspired by BIDS)

[valosekj]

On the left the t1 image, on the right the t2 image

please avoid using t1 t2 for times points 1-2, because in MRI we often use t1 and t2 for the contrasts. Instead, you could use ses-1, ses-2 (inspired by BIDS)

Or maybe, ses-M0 and ses-M12 could be used to be consistent with git-annex/canproco.

[plbenveniste]

As suggested in this comment, I tried using sct_label_utils -remove-reference to get a common list of vertebral levels in both segmentation file, in order to use them during registration. I apply it to both images using the other as a reference. I also overwrite the images.

os.system('sct_label_utils -i ' + first_vert_output_path + ' -remove-reference ' + second_vert_output_path + ' -o ' + first_vert_output_path)
t1 = time()
print("time to remove-reference",t1-t0)
os.system('sct_label_utils -i ' + second_vert_output_path + ' -remove-reference ' + first_vert_output_path + ' -o ' + second_vert_output_path)
t2 = time()
print("time to remove-reference",t2-t1)

However, I am surprised by a very long computation time for the images (which is much longer than the previous method which consisted of loading the data, removing outliers and saving it.

For the first image: 147 seconds. For the second image: 161 seconds

Terminal output:

[jcohenadad]

about #38 (comment): it is surprising indeed, given the triviality of the task. Could you please open an issue on the SCT repos? Thanks!