dcm2mnc assigns wrong Slice indices to 4D PET dicom series

[jpcoutu]

I have encountered a 4D PET DICOM series that gets reconstructed properly using an old dcm2mnc (2014, minc-2.2.00, unsure which repo/commit), and not using the current develop branch. I think this is due to a problem within get_file_indices() in dicom_read.c.

Running dcm2mnc in debug mode, I can see that new indices are being added to the Slice list for the same coordinate value, while they should only be added to the Time list (which they also are). I have a total of 47 slices (therefore max index of 47), and I run out of new Slice indices a quarter of the way through (since I have 4 time points). I have pasted the relevant section of the dcm2mnc standard output below and bolded the problematic part.

This seems due to the fact that the Instance Number goes up first via the time axis rather than the slice axis (I have pasted a few dcmdump fields below the debug output to show this). Seems maybe counter intuitive in terms of acquisition, but should be able to reconstruct this nevertheless, as the old version of dcm2mnc did.

I am wondering whether that old version was assigning slice indices by coordinate values (similar to the way it still seems to be done for IMA file types)? The G.file_type I get (as parsed by dcm2mnc) for these files is N4DCM. I am happy to attempt fixing this issue, but would appreciate any feedback on the best course of action to do so. Please let me know if any other information would be useful to share.

Debug output:

File fail/171
0. Slice axis length 47

1. Echo axis length 1
   (0020,0105)=-1
   (0054,0101)=4
2. Time axis length 4
3. Phase axis length 1
4. ChmSh axis length 1
   get_coordinate_info
   dircos 0.000000 -0.000000 1.000000 -0.000000 -1.000000 -0.000000 -1.000000 -0.000000 0.000000
   Volume_to_world slice=Z row=Y column=X
   Orientation is TRANSVERSE
   Using (0018,0050) for slice thickness
   Swapping direction of Row Y
   Swapping direction of Column X
   coordinate 127.000000 127.000000 -150.399994, start 127.000000 127.000000 -150.399994
5. Z axis length 47 step 3.270, start -150.400, cosines 0.000,-0.000,1.000
6. Y axis length 128 step -2.000, start 127.000, cosines 0.000,1.000,0.000
7. X axis length 128 step -2.000, start 127.000, cosines 1.000,0.000,-0.000
   SOP Class UID: 1.2.840.10008.5.1.4.1.1.128
   Images in acquisition: -1
   Acquisitions in series: -1
   3D raw partitions: -1

Using GE-specific PET information.
Pixel minimum -32768.0000000000 maximum 32767.0000000000
Window minimum -20521.9102720000 maximum 20521.2839930000
First Slice at 1,-150.399994,0.000000
First Echo at 1,0.000000,0.000000
First Time at 1,0.000000,300.000000
First Phase at 1,0.000000,0.000000
First ChmSh at 1,0.000000,0.000000

File fail/126
get_coordinate_info
dircos 0.000000 -0.000000 1.000000 -0.000000 -1.000000 -0.000000 -1.000000 -0.000000 0.000000
Volume_to_world slice=Z row=Y column=X
Orientation is TRANSVERSE
Using (0018,0050) for slice thickness
Swapping direction of Row Y
Swapping direction of Column X
coordinate 127.000000 127.000000 -150.399994, start 127.000000 127.000000 -150.399994
Need to add index 2 to Slice list, 1/47
Added Slice coordinate -150.399994 at 2 127.000000 127.000000 -150.399994
Need to add index 2 to Time list, 1/4
Added Time coordinate 300.000000 at 2 127.000000 127.000000 -150.399994

File fail/109
get_coordinate_info
dircos 0.000000 -0.000000 1.000000 -0.000000 -1.000000 -0.000000 -1.000000 -0.000000 0.000000
Volume_to_world slice=Z row=Y column=X
Orientation is TRANSVERSE
Using (0018,0050) for slice thickness
Swapping direction of Row Y
Swapping direction of Column X
coordinate 127.000000 127.000000 -150.399994, start 127.000000 127.000000 -150.399994
Need to add index 3 to Slice list, 2/47
Added Slice coordinate -150.399994 at 3 127.000000 127.000000 -150.399994
Need to add index 3 to Time list, 2/4
Added Time coordinate 600.000000 at 3 127.000000 127.000000 -150.399994

File fail/123
get_coordinate_info
dircos 0.000000 -0.000000 1.000000 -0.000000 -1.000000 -0.000000 -1.000000 -0.000000 0.000000
Volume_to_world slice=Z row=Y column=X
Orientation is TRANSVERSE
Using (0018,0050) for slice thickness
Swapping direction of Row Y
Swapping direction of Column X
coordinate 127.000000 127.000000 -150.399994, start 127.000000 127.000000 -150.399994
Need to add index 4 to Slice list, 3/47
Added Slice coordinate -150.399994 at 4 127.000000 127.000000 -150.399994
Need to add index 4 to Time list, 3/4
Added Time coordinate 900.000000 at 4 127.000000 127.000000 -150.399994

File fail/159
get_coordinate_info
dircos 0.000000 -0.000000 1.000000 -0.000000 -1.000000 -0.000000 -1.000000 -0.000000 0.000000
Volume_to_world slice=Z row=Y column=X
Orientation is TRANSVERSE
Using (0018,0050) for slice thickness
Swapping direction of Row Y
Swapping direction of Column X
coordinate 127.000000 127.000000 -147.129990, start 127.000000 127.000000 -147.129990
Need to add index 5 to Slice list, 4/47
Added Slice coordinate -147.129990 at 5 127.000000 127.000000 -147.129990

File fail/176
get_coordinate_info
dircos 0.000000 -0.000000 1.000000 -0.000000 -1.000000 -0.000000 -1.000000 -0.000000 0.000000
Volume_to_world slice=Z row=Y column=X
Orientation is TRANSVERSE
Using (0018,0050) for slice thickness
Swapping direction of Row Y
Swapping direction of Column X
coordinate 127.000000 127.000000 -147.129990, start 127.000000 127.000000 -147.129990
Need to add index 6 to Slice list, 5/47
Added Slice coordinate -147.129990 at 6 127.000000 127.000000 -147.129990

dcmdump values for those 6 files above (total of 47 slices, 4 time points)

dcmdump +P 0020,0013 +P 0054,1330 +P 0020,1041 fail/171
(0020,0013) IS [1] # 2, 1 InstanceNumber
(0054,1330) US 1 # 2, 1 ImageIndex
(0020,1041) DS [-150.39999389648] # 16, 1 SliceLocation
dcmdump +P 0020,0013 +P 0054,1330 +P 0020,1041 fail/126
(0020,0013) IS [2] # 2, 1 InstanceNumber
(0054,1330) US 48 # 2, 1 ImageIndex
(0020,1041) DS [-150.39999389648] # 16, 1 SliceLocation
dcmdump +P 0020,0013 +P 0054,1330 +P 0020,1041 fail/109
(0020,0013) IS [3] # 2, 1 InstanceNumber
(0054,1330) US 95 # 2, 1 ImageIndex
(0020,1041) DS [-150.39999389648] # 16, 1 SliceLocation
dcmdump +P 0020,0013 +P 0054,1330 +P 0020,1041 fail/123
(0020,0013) IS [4] # 2, 1 InstanceNumber
(0054,1330) US 142 # 2, 1 ImageIndex
(0020,1041) DS [-150.39999389648] # 16, 1 SliceLocation
dcmdump +P 0020,0013 +P 0054,1330 +P 0020,1041 fail/159
(0020,0013) IS [5] # 2, 1 InstanceNumber
(0054,1330) US 2 # 2, 1 ImageIndex
(0020,1041) DS [-147.12998962402] # 16, 1 SliceLocation
dcmdump +P 0020,0013 +P 0054,1330 +P 0020,1041 fail/176
(0020,0013) IS [6] # 2, 1 InstanceNumber
(0054,1330) US 49 # 2, 1 ImageIndex
(0020,1041) DS [-147.12998962402] # 16, 1 SliceLocation

[vfonov]

I hope @rdvincent can help with this

[jpcoutu]

I have created a temporary fix to this problem. I say temporary because I don't think this is the right way to go about it, but it does what needs to be done without creating a problem for other dicom series (following our testing). Basically this is a sub called fix_dimensions() that is called right before sort_dimensions() in dicom_to_minc.c. It fixes the problem mentioned above, and another one that I have not posted (different dicom series), which is the same issue but for indices of the Time list getting the wrong assignment because of a similar situation. Both problematic cases were 4D PET GEMS dicom series.

static void
fix_dimensions(File_Info *fi_ptr, General_Info *gi_ptr)
{
    Mri_Index imri;
    int found, duplicate, unique;
    int i, j, k;

    if (gi_ptr->max_size[SLICE] > 1 && gi_ptr->max_size[TIME] > 1) {
        for (imri = 0; imri < MRI_NDIMS; imri++) {
            if (gi_ptr->cur_size[imri] > 1 && !((G.file_type == N4DCM)
                && (imri == TIME) && (gi_ptr->num_slices_in_file > 1))) {
                duplicate = 0;
                unique = 1;
                for (i = 0; i < gi_ptr->max_size[imri] - 1; i++) {
                    for (j = i + 1; j < gi_ptr->max_size[imri]; j++) {
                        if (gi_ptr->coordinates[imri][i] == gi_ptr->coordinates[imri][j]) {
                            duplicate = 1;
                        }
                        else {
                            unique = 0;
                        }
                    }
                }
                /* If we found a duplicate in this context, the fix is needed,
                   unless all files have the same coordinate */
                if (duplicate == 1 && unique == 0) {
                    k = 0;
                    for (i = 0; i < gi_ptr->num_files; i++) {
                        found = 0;
                        for (j = 0; j < k; j++) {
                            if (fi_ptr[i].coordinate[imri] == gi_ptr->coordinates[imri][j]) {
                                found = 1;
                            }
                        }
                        if (found == 0) {
                            k++;
                            gi_ptr->coordinates[imri][k-1] = fi_ptr[i].coordinate[imri];
                        }
                        fi_ptr[i].index[imri] = k;
                    }
                }
            }
        }
    }
}

[jpcoutu]

Since my last post I encountered two other variants of this case, both GEMS PET, and I figured I would share them and post an updated fix_dimensions() sub. The first case had no obvious field to determine the number of time frames, so it was needed to allow the time dimension to grow (it says close to line 900 in dicom_read.c that time dimension is allowed to grow, but isn't). I enabled this for some cases:

          if (imri != SLICE && gi_ptr->max_size[imri] <= 1 && !(imri == TIME &&
                !((G.file_type == N4DCM) && (gi_ptr->num_slices_in_file > 1)))) {
            <...>
            continue;
          }

This however created an issue where the second dicom file processed through get_file_indices() would not get range-checked since max_size[TIME] had not yet grown above 1. Added a fix below in the fix_dimensions sub.

The other case had a bogus dicom field InstanceNumber, which took precedence over the field ImageIndex in get_file_indices(). When I gave ImageIndex precedence, it provided the appropriate reconstruction. My previous fix_dimensions caught that, but fixed it wrong (likely because of bad initial sorting with the InstanceNumber). Below is the updated sub fixing both these new cases. Happy to pull-request this, but I think maybe something different needs to happen in dcm2mnc specifically for GEMS PET to fix these issues. Maybe ACR_PET_Image_index (ImageIndex) could override ACR_Temporal_position_identifier for these cases for instance (but that may just fix one of the many issues). Happy to discuss further by email or in person.

static void
fix_dimensions(File_Info *fi_ptr, General_Info *gi_ptr)
{
    Mri_Index imri;
    int found, duplicate, unique, count, expected;
    int i, j, k;

    /* Hopefully temporary fix to bad index assignment for 4D PET volumes,
       sometimes for the slice coordinate, sometimes for the time coordinate  */
    if (gi_ptr->max_size[SLICE] > 1 && gi_ptr->max_size[TIME] > 1) {
        for (imri = 0; imri < MRI_NDIMS; imri++) {
            if (gi_ptr->cur_size[imri] > 1 && !((G.file_type == N4DCM)
                && (imri == TIME) && (gi_ptr->num_slices_in_file > 1))) {
                duplicate = 0;
                unique = 1;
                for (i = 0; i < gi_ptr->max_size[imri] - 1; i++) {
                    for (j = i + 1; j < gi_ptr->max_size[imri]; j++) {
                        if (gi_ptr->coordinates[imri][i] == gi_ptr->coordinates[imri][j]) {
                            duplicate = 1;
                        }
                        else {
                            unique = 0;
                        }
                    }
                }
                /* If we found a duplicate in this context, the fix is needed,
                   unless all files have the same coordinate */
                if (duplicate == 1 && unique == 0) {
                    k = 0;
                    for (i = 0; i < gi_ptr->num_files; i++) {
                        found = 0;
                        for (j = 0; j < k; j++) {
                            if (fi_ptr[i].coordinate[imri] == gi_ptr->coordinates[imri][j]) {
                                fi_ptr[i].index[imri] = j+1;
                                found = 1;
                                break;
                            }
                        }
                        if (found == 0) {
                            k++;
                            gi_ptr->coordinates[imri][k-1] = fi_ptr[i].coordinate[imri];
                            fi_ptr[i].index[imri] = k;
                        }
                    }
                }
            }
        }
    }
    /* A consequence of allowing the time dimension to grow is that sometimes the
       second file in a series will not be range-checked because max_size hasn't
       grown above one yet, creating a second index that shouldn't be. Removing the
       max_size[TIME] > 1 condition on the range check doesn't work in cases where
       the fix above is needed, therefore this following fix is needed. */
    if (gi_ptr->max_size[TIME] == 2) {
        count = 0;
        expected = 0;
        for (i = 0; i < gi_ptr->num_files; i++) {
            if (fi_ptr[i].index[TIME] == gi_ptr->indices[TIME][1]) {
                count++;
            }
            else {
                expected++;
            }
        }
        /* If we found only one file with a different index, fix is needed */
        if (count == 1 && expected > 1) {
            gi_ptr->cur_size[TIME] = 1;
            gi_ptr->max_size[TIME] = 1;
            for (i = 0; i < gi_ptr->num_files; i++) {
                fi_ptr[i].index[TIME] = gi_ptr->indices[TIME][0];
            }
        }
    }
}

[jpcoutu]

This is all fixed/included in #113