vfonov
commented
6 years ago I hope @rdvincent can help with this
Open jpcoutu opened 6 years ago
vfonov
commented
6 years ago I hope @rdvincent can help with this
jpcoutu
commented
6 years ago 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;
}
}
}
}
}
}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];
}
}
}
}This is all fixed/included in https://github.com/BIC-MNI/minc-tools/pull/113
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
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 SliceLocationdcmdump +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 SliceLocationdcmdump +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 SliceLocationdcmdump +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 SliceLocationdcmdump +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 SliceLocationdcmdump +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