Mirror issue

[SitronX]

Hello @mlavik1,
I have encountered very particular issue. When dataset comes from CT machine, i guess it can be scanned two ways. The patient is either scanned from head to feet, or from feet to head.

The major issue is, when the patient is scanned from feet to head, the library doesnt know that and after careful observation, all organs are mirror flipped (organs that should be in left part of the body are in right part).

I also havent noticed that before, but surgeons noticed :)

The simple temporary fix for these specific datasets is to reverse the data texture array, before the texture is generated. Like this (same for gradient texture)

It will then correctly show the dataset, maybe rotated around little bit, but organs are not mirror flipped anymore and it matches.

The issue is, this is needed only for CT scans that are from feet to head, CT scans from head to feet are probably rendering correctly (I dont have head to feet datasets).

The library should ideally detect how are the data stored and then subsequently flip the array if needed. I guess comparing two slice positions and recognizing which type of scan it is should work.

Here is how the dataset looks in Slicer3D for comparison

FeetToHead.mp4

This is how it looks in this version of the library, organs are basically mirror flipped

This is after reversing the texture array, organs are now in correct place similar to Slicer3D

The pictures are only for position comparisons, i probably have slightly different density interval on both pictures.

[SitronX]

So i have been searching for couple of hours and i have a problem to really find any DICOM dataset that starts from top-to-bottom. I only found one dataset, which was behaving really odly with file ordering, so it might have not been even top-to-bottom and it was sorted artificially. For simpleITK image sequence importer this should work.

IsHeadFeetDataset method is here

private bool IsHeadFeetDataset(Image firstImage,Image lastImage)
{
    if (TryGetPositionInPatient(firstImage, out Vector3 firstPosition))
    {
        if (TryGetPositionInPatient(lastImage, out Vector3 secondPosition))
        {
            if (firstPosition.z > secondPosition.z)                  
                return true;                   
        }
    }
    return false;
}

TryGetPositionInPatient method is here

private bool TryGetPositionInPatient(Image sliceImage,out Vector3 position)
{
    try
    {
        List<string> metadataKeys = sliceImage.GetMetaDataKeys().ToList();

        if (metadataKeys.Contains("0020|0032"))
        {
            string imagePositionPatient = sliceImage.GetMetaData("0020|0032");    //0020|0032 tag for getting location
            string[] arr = imagePositionPatient.Split('\\');

            float x = float.Parse(arr[0], CultureInfo.InvariantCulture);
            float y = float.Parse(arr[1], CultureInfo.InvariantCulture);
            float z = float.Parse(arr[2], CultureInfo.InvariantCulture);

            position = new Vector3(x, y, z);
            return true;
        }
        position = Vector3.zero;
        return false;
    }
    catch
    {
        position = Vector3.zero;
        return false;
    }
}

For simpleITK file importer it should be something like this. The issue is, i cannot really read slice position metadata from nrrd.

This is a showcase for SimpleITK importer, something like this should be probably added to every importer if it is having same mirror problems.

[mlavik1]

Oh, wow that's actually a pretty big issue. I hadn't noticed it all this time. Thanks for reporting!

And yes, I think that would be a good solution!
Should I make a PR, or do you want to create one from your fork?
I can add a similar change to OpenDICOM as well.

[SitronX]

@mlavik1
ImageFileSequenceImporter should work overall, i dont have proper head-to-feet datasets, but the logic should be correct.

After further testing, the ExtractSlice method does work, it extracts the slice but i cant really get it to extract any metadata from the specific slice. The metadata array is always empty. The ImageFileReader contains overall metadata, but that is not really helpful. Do you know the way how to extract the specific slice metadata from formats like nrrd or nifti? Is it even possible?

Also, it is be possible that sometimes the location tag might be missing, I think it is possible even in DICOM sequence. This is a real pain, because we would probably have to scan all metadata and hope there is some clue. I wonder how other software like lets say Slicer3D really do it.
Edit: I was wrong, the 0020|0032 is required tag in DICOM

My temporary suggestion: If it contains the position metadata and the flip is required, lets flip it. If there is no metadata, it will not do anything with possibility that it might be incorrectly mirrored. I think majority of the datasets are from feet-to-head, so lets default to that. When somebody loads the head-to-feet dataset with no location metadata, it will mirror, we can add option to flip it in realtime to manually correct it.

This would be temporary solution. We should probably still find better way to do it, if it is possible.

Well I could probably make PR if the metadata extraction from files like nrrd an nifti worked. So it is still an open issue.

[mlavik1]

Hi @SitronX

Yes, I agree with your suggestion.
And I think you're right that the majority of datasets are feet-to-head, because according to this it would seem like that is the default: "The z-axis is increasing toward the head of the patient.".

I'm not sure what we could do if the position tag is missing.. My first thought was to use the location tag, but apparently that is " relative to an unspecified implementation specific reference point". I guess if they scanned a patient top-to-bottom without adding a position element then that is probably a mistake, that would cause issues in other software as well? If we find a head-to-bottom dataset we could probably verify that by removing the position tag and trying to import it elsewhere.
But for now we can probably just assume that. Anything's better than the current inverted bodies 😅

I have a bunch of DICOM datasets I can test with as well, and also some NRRD/NIFTI. I'll see if I can have a look at it tomorrow or in the weekend. After that I'll be away for a week :)

[SitronX]

@mlavik1
Hello, do you got any update on this topic? Have you managed to somehow read metadata from NRRD in your tests?

[mlavik1]

Hi @SitronX !
I had a busy first week back at work after my holiday, so I haven't had any time to work on this project, but I'm planning to do a bit today and more on Tuesday/Wednesday. I'll let you know if I get somewere with it :)

I'll also start creating a PR based on the changes you posted above (saw you pushed them to your repo as well).

[SitronX]

Hello @mlavik1
Since you looking to my code above. Just a quick correction as i stumbled upon it back then.

The sequenceSeries reverse as i showed above in SimpleITKImageSequenceImporter.cs is wrong. It needs to be read normally and then the whole pixelData array at the end must be reversed similarly as it is done in SimpleITKFileImporter.cs. Otherwise the data will be little bit corrupted due to wrong pixel ordering between the slices.

[mlavik1]

Ok, thanks for the info! :)

[mlavik1]

Hi again!
So I've started preparing a branch with the changes you suggested, and done some testing and reading.
So I have one question: According to this the Z axis should always be increasing towards the head.
And according to this DICOM should use LPS (Left, Posterior, Superior) coordinate system.
So should it really be necessary to check the positions of each slice? Or can we just assume that all DICOM datasets are in LPS system with Z values increasing towards the head?
Then we just need to convert from LPS to Unity's coordinate system (x=-x, y=z, z=-y), or just reverse the pixel array as you did (which I think should have the same effect, except maybe the rotation will be different?)

And regarding NRRD: It seems to have a "space" field: https://teem.sourceforge.net/nrrd/format.html. Could that be what we're looking for?
I suppose SimpleITK might already handle this for us though..
To get the same orientation for NRRD and DICOM, maybe we could call SimpleITK.DICOMOrient(image, "LPS")?

Do you have some NRRD files that have a different orientation / coordinate system, or are all your NRRD datasets similar as well?

I pushed my suggested changes here: https://github.com/mlavik1/UnityVolumeRendering/tree/coordinate-systems
(still WIP).
Does that work for your NRRD/NIFTI datasets as well? (I need to find some better ones where it's easier to see if they are mirrored or not :D)

[SitronX]

Hello @mlavik1

According to this the Z axis should always be increasing towards the head.

Yes i think that is right. With Z, head should always have higher vale than feet. I havent found anything saying otherwise.

And according to this DICOM should use LPS (Left, Posterior, Superior) coordinate system.

Aha, you are right. Supposedly, DICOM should always be stored in LPS. If you view DICOM in Slicer3D, Slicer3D specifically
converts it to RAS for its own internal usage (i guess). But overall, dicom on disc should be always in LPS, that is good to know. Althought i dont think the RAS would cause the issue, cause the Z also increases towards the head of the patient.

So should it really be necessary to check the positions of each slice?

In my code above, not every slice is checked for position. Only first and last slice of dataset is checked and Z is compared, to determine the type of acquisition.

While the LPS should always be the case, i dont think it actually tells the acquisition order. I think it is just a coordinate system, and both head-to-feet or feet-to-head dataset can use same coordinate system, while acquisition order would be different. It would still mean, that head has higher Z value than feet, thus it is increasing towards the head.

But without checking Z value, either head-to-feet or feet-to-head would still cause mirror issue (depending if it is flipped defaulty).

Here are described possible types of scans,

I have also found this interesting post from user, where he was dealing with same problem, this was sent to him. Which describes to him how it should be dealt with. The lady mentiones there, that these tags are the only reliable way how to represent dataset correctly (in right acquisition order). It is really interesting read, i recommend it.

It is also possible to get acquisition order from the already mentioned tag, but this tag is not required and might be missing sometimes. Reading Z value from slice position in patient is reliable, because the position tag should always be there.

or just reverse the pixel array as you did (which I think should have the same effect, except maybe the rotation will be different?)

Yes you are right, the orientation is slightly different and i had to correct it. Switching axis like you said (x=-x, y=z, z=-y) would probably be better, if it works :D

So i have checked your solution, and you have it backwards, cause my nrrd mirrors :D

You have
bool isDatasetReversed=SimpleItkUtils.IsHeadFeetDataset(firstSlice,lastSlice);

and it should be

bool isDatasetReversed=!SimpleItkUtils.IsHeadFeetDataset(firstSlice,lastSlice);

The data reverse should happen, if it is feet-to-head dataset. (Edit: Just noticed you have it switched only in file reader, sequence reader is alright)

But i havent checked the rest of the code, i only tested the one nrrd i have :D
For testing, the ircad datasets here are nice. It is DICOM only and not nrrd, but it can be converted in Slicer3D for your tests. They are really detailed and nice datasets.

[mlavik1]

Hi @SitronX

Thanks for the extra info! That mailing list reply you linked was really informative :)

I think it is just a coordinate system, and both head-to-feet or feet-to-head dataset can use same coordinate system, while acquisition order would be different.

Ah, yes, acquisition order might be in any direction. So if we simply used the SimpleItk DICOM reader to read the dataset slice by slice and then manually tried to construct a dataset (like that person asking in the mailing list), we would need to order the slices ourselves. However, since we're using the ImageSeriesReader class (which constructs a 3D pixel array from a series of DICOM files) I would assume that it uses the position tag (if present) to ensure that the resulting 3D pixel array will have the Z value increase towards the head (bottom-up)? Or am I still misunderstanding something? 😅 (I have a feeling that I am haha)

So i have checked your solution, and you have it backwards, cause my nrrd mirrors :D

Ooops! You're right haha. I'l lfix that right away :D

But i havent checked the rest of the code, i only tested the one nrrd i have :D
For testing, the ircad datasets here are nice. It is DICOM only and not nrrd, but it can be converted in Slicer3D for your tests. They are really detailed and nice datasets.

Thanks, I'll try out those datasets!
And I didn't know Slicer3D could convert to NRRD. I guess that might be a good way to test, so I can test the "same" datasets through DICOM and NRRD 😁

[SitronX]

haha sorry for the confusion :D

You are right. So supposedly the GetGDCMSeriesFileNames() handles the ordering. It is not exactly specified how, but it says it orders it with several strategies. So that means, that the ordering should be done automatically for us with SimpleITK and ImageSeriesReader.

Since SimpleITK handles this ordering in ImageSeriesReader, is it possible that something similar happens with FileImporter, where it loads it automatically in same way (feet first), no matter what type of scan it is? So that would mean the only required thing is to defaultly flip it and it should always be correct no matter the scan type, cause SimpleITK should always order it as feet first.

But hey, atleast here is a path how to fix the standard dicom loader, if SimpleITK is not used :D

Again, sorry for the confusion. I just wish the itk/simpleITK documentation was more detailed on this topic.

[mlavik1]

Hi again @SitronX !

Oh, right! That might be the function that does it then :)

Since SimpleITK handles this ordering in ImageSeriesReader, is it possible that something similar happens with FileImporter, where it loads it automatically in same way (feet first), no matter what type of scan it is?

Yes, I think that's the case!
I did a quick experiment with an NRRD file. I opened (with a hex viewer - because a text viewer can end up breaking it) and changed "space: left-posterior-superior" to "space: left-anterior-superior". When opening them in 3D Slicer the original and modified files look different, which is expected since we changed the coordinate space without changing the data. However, in the Unity volume rendering plugin nothing changed - unless I call SimpleITK.DICOMOrient(image, "LPS").

So it might seem like SimpleITK reads the coordinate space, but only converts to LPS on demand, so we should probably do that :)

But hey, atleast here is a path how to fix the standard dicom loader, if SimpleITK is not used :D

Yes, that's right! It already reads the position tag and uses that to order the slices (see CalcSliceLocFromPos) but the way it's done is not right... the slices will be sorted correctly relative to each other, but whether the order is reversed or not will "random" 😅 Thanks to you we know how to fix that now :D

I just wish the itk/simpleITK documentation was more detailed on this topic.

Yes, it's hard to find out exactly what's going on.

So the solution to this issue would then be to:

• Execute SimpleITK.DICOMOrient(image, "LPS")
• Always reverse the pixel array, or convert the coordinate system some other way.
• Fix the OpenDICOM importer's slice ordering

Is that right?

Regarding the second one, reversing the pixel array is definitely the simplest fix. We might want to swap it in place though, instead of how it's done now, to avoid extra memory allocations (these arrays can be big! haha)
Alternatively we could just import the dataset as before, and instead change the scale of the mesh container GameObject. Not sure what's the nicest solution to be honest 😅

[SitronX]

So the solution to this issue would then be to:

Yep, that sound great.

Alternatively we could just import the dataset as before, and instead change the scale of the mesh container GameObject.

Ahh, that is actually great solution. Changing the x scale to -x of the rendered object. This is so much better than reversing the array and so simple. I am gonna ditch reversing the array in my project and use this instead :D

[mlavik1]

Ahh, that is actually great solution. Changing the x scale to -x of the rendered object. This is so much better than reversing the array and so simple. I am gonna ditch reversing the array in my project and use this instead :D

Before you celebrate: Doing that will apparently break the slicing plane tool 😅
I think it's because the slicing plane is relative to the outer volume GameObject - we should probably make it relative to the "VolumeContainer" object instead! The actual plane can still be attached to the outer GameObject, but we need to change GetMatrix in CrossSectionPlane.cs I think..

I'll make a new branch with the changes we talked about, so we can eventually create a PR :)

[mlavik1]

Ok, I think the PR should be ready for review now: #170

It has a lot of changes, so if there's something there you disagree with or think should be different please let me know 😁
And no hurry!

I think the SimpleITK importer is the most important to fix, and then I can do a separate PR for fixing the OpenDICOM importer.

[mlavik1]

Fixed for the Open DICOM importer as well!