Normalizing a group of images for the same patient

[SamuelJoutard]

Hi,

I was wondering how intensity normalizations work when a patient contains multiple images. For instance when we do ZNormalization, are the normalization statistics (i.e. means and std) computing jointly on all images or per image.

To give some more details I see two different use case example which would require two different settings. First, if the images of the patients correspond to different modalities then you obviously want to compute those statistics per image. On the other hand, when you do unimodal longitudinal registration, you want intensities to be normalized jointly in order to be able to match intensities.

Thank you,

Samuel

[ReubenDo]

Hi Samuel,

The current normalization is per image.

For unimodal registration, wouldn't it still sometimes make more sense to normalize each image independently? For example, if there is an intensity shift between the two images, only an image-specific normalization would work, no?

Reuben

[fepegar]

Hi all,

I agree with @ReubenDo, I was just plotting some histograms to show what happens when you normalize separately or jointly.

[fepegar]

Initial values:

After ZNorm (top: independent stats; bottom: joint stats):

(Note the second mode is not centered on 0 because I didn't use foreground masks to compute the values)

[fepegar]

Here's a version using values above the mean (i.e. foreground) to compute the stats (top: independent stats; bottom: joint stats):

[fepegar]

And these are the images I'm using (top: blue, bottom: orange):

[SamuelJoutard]

Hi all,

Indeed, I agree with you Reuben in case of intensity shift it should be done per image at least for the mean part.

The thing is if you use SSD as similarity measure, and if images are acquired with the same scanner then you definitely want to do joint normalization because otherwise unchanged organs between the two images will have different intensities.

Thank you for those illustrations Fernando. I still think that there are some cases where you want to compute the stats jointly especially when the key information for your loss or your processing is the relationship between the intensities of the two images (or the registered images).

[SamuelJoutard]

Maybe this is too specific to some use cases so I did not mean to suggest it as a feature request but I guess it could be specified in the documentation in case some people expect joint normalization.

[fepegar]

I still think that there are some cases where you want to compute the stats jointly especially when the key information for your loss or your processing is the relationship between the intensities of the two images (or the registered images).

Ok, then this looks like a feature request :)

What do you think would be a good API for this? Would something like this work?

transform = torchio.ZNormalization(joint_stats=True)

[fepegar]

Maybe this is too specific to some use cases so I did not mean to suggest it as a feature request but I guess it could be specified in the documentation in case some people expect joint normalization.

It is specific, but I think it makes sense as a feature to add. I'm biased towards (brain) segmentation, so it's good to get input from people working on different tasks.

[SamuelJoutard]

Alright cool then what you suggested looks like a good API for me.

If you want to handle all use cases then joint_stats could also be a dictionary with groups of images to be jointly normalized per subjects. For instance for longitudinal registration using multiple modalities, if your subject is:

subjects_list = [
    torchio.Subject(
        moving_MRI=torchio.Image(path_MRI[i], torchio.INTENSITY),
        moving_CT=torchio.Image(path_CT[i], torchio.INTENSITY),
        label_m=torchio.Image(path_lbl[i], torchio.LABEL),
        reference_MRI=torchio.Image(path_MRI[j], torchio.INTENSITY),
        reference_CT=torchio.Image(path_CT[j], torchio.INTENSITY)
        label_r=torchio.Image(path_lbl[j], torchio.LABEL),
    )
    for i in range(N_train)
    for j in range(i+1, N_train)
]

Then the joint_stats could accept a boolean (to either do per image normalization or to jointly normalize all patient images) or a dict to be a bit more modular like:

joints_stats = {
    0: ["moving_MRI", "reference_MRI"],
    1: ["moving_CT", "reference_CT"],
}

Where the normalization would be perform jointly per group.

That being said, the API you suggested is all good for me so this is just a suggestion.

Thank you for your help.

[fepegar]

(@SamuelJoutard I edited your comment to improve a bit the readability)

Yeah, that makes sense. I'll think of a way to do this and hopefully work something out by Friday.

[SamuelJoutard]

Thanks a lot for your help!

[fepegar]

I'll leave this open, as a reminder.

[fepegar]

Apologies for the delay, it's not trivial to solve this issue.

[fepegar]

If a) your images have the same shape and b) you haven't given up on using TorchIO, I think you can now make this work, as we now support 4D images and at the moment the stats are computed from the whole image. Let me know if you want to follow that path.

[SamuelJoutard]

Hi Fernando,

Thanks you for not giving up on this and proposing this solution. At the moment, my images do not have the same shape but I could consider some preprocessing to standardize this. I am of course still using TorchIO for the many others convenient functionality proposed in your library.

I will let you know if I end up using this new functionality and provide feedbacks if possible. Thank you for the help!

Samuel

[fepegar]

Ok, great! You can use the Resample transform and pass a reference image or an affine transform.

[fepegar]

Closing this for now. We can reopen if there's still interest.