changeset 2260:e744f0414183

Improve support for DTI
author bert <bert>
date Sun, 09 Apr 2006 15:39:04 +0000
parents 915eb0cdc0df
children 8a8ece4eeb00
files conversion/dcm2mnc/dicom_to_minc.c
diffstat 1 files changed, 267 insertions(+), 3 deletions(-) [+]
line wrap: on
line diff
--- a/conversion/dcm2mnc/dicom_to_minc.c
+++ b/conversion/dcm2mnc/dicom_to_minc.c
@@ -8,7 +8,10 @@
    @CREATED    : January 28, 1997 (Peter Neelin)
    @MODIFIED   : 
    * $Log: dicom_to_minc.c,v $
-   * Revision 1.17  2006-02-09 20:54:29  bert
+   * Revision 1.18  2006-04-09 15:39:04  bert
+   * Improve support for DTI
+   *
+   * Revision 1.17  2006/02/09 20:54:29  bert
    * More changes to dcm2mnc
    *
    * Revision 1.16  2005/12/13 17:31:13  bert
@@ -174,7 +177,7 @@
    provided "as is" without express or implied warranty.
    ---------------------------------------------------------------------------- */
 
-static const char rcsid[] = "$Header: /private-cvsroot/minc/conversion/dcm2mnc/dicom_to_minc.c,v 1.17 2006-02-09 20:54:29 bert Exp $";
+static const char rcsid[] = "$Header: /private-cvsroot/minc/conversion/dcm2mnc/dicom_to_minc.c,v 1.18 2006-04-09 15:39:04 bert Exp $";
 #include "dcm2mnc.h"
 
 const char *World_Names[WORLD_NDIMS] = { "X", "Y", "Z" };
@@ -676,6 +679,21 @@
     if (strstr(str_ptr, "MOSAIC") != NULL)
         return 1;               /* slam dunk, this is mosaic data */
 
+    /* For some reason, some Localizer scans look like MOSAICs in that
+     * they appear to be upsampled in such a way that their
+     * rows/columns are larger than the acquisition matrix. The
+     * details of this are probably buried in the Siemens ASCCONV dump
+     * and therefore I have not been able to figure out how to tweak
+     * the MOSAIC stuff to get this right. Instead I am taking the
+     * cheesy way out and declaring that anything that is NOT
+     * explictly marked MOSAIC, and IS marked a localizer, should not
+     * be treated as a mosaic.
+     */
+    str_ptr = acr_find_string(group_list, ACR_Series_description, "");
+    if (strstr(str_ptr, "localizer") || strstr(str_ptr, "LOCALIZER")) {
+        return 0;
+    }
+
     /* OK, we did not find the word "mosaic" in the image type.  But this
      * could still be a mosaic image.  Let's look for some other clues.
      */
@@ -688,6 +706,114 @@
     return 0;                   /* probably not mosaic */
 }
 
+#define MAXVM 32
+
+Acr_Group
+parse_siemens_proto2(Acr_Group group_list, Acr_Element element)
+{
+    int byte_cnt;
+    int byte_pos;
+    unsigned char *byte_ptr;
+    int n_items;
+    int n_values;
+    int vm;
+    char vr[4];
+    int syngodt;
+    int i;
+    char name[64];
+    char *value[MAXVM];
+    int len;
+
+    byte_cnt = element->data_length;
+    byte_ptr = element->data_pointer;
+    byte_pos = 0;
+
+    /* See if the magic 8-byte header is present. If not, we start
+     * right off with the number of elements.
+     */
+    if (byte_ptr[0] == 'S' && byte_ptr[1] == 'V' &&
+        byte_ptr[2] == '1' && byte_ptr[3] == '0') {
+        byte_pos += 8;          /* Skip header */
+    }
+
+    n_items = *((int *)(byte_ptr + byte_pos));
+
+    byte_pos += 8;              /* Skip # items and 4 bytes of unknown junk */
+
+    while (n_items-- > 0) {
+        strncpy(name, (byte_ptr + byte_pos), 64);
+        byte_pos += 64;
+
+        vm = *((int *)(byte_ptr + byte_pos));
+        byte_pos += 4;
+
+        strncpy(vr, (byte_ptr + byte_pos), 4);
+        byte_pos += 4;
+
+        syngodt = *((int *)(byte_ptr + byte_pos));
+        byte_pos += 4;
+
+        n_values = *((int *)(byte_ptr + byte_pos));
+        byte_pos += 4;
+
+        byte_pos += 4;          /* skip dummy */
+
+        if (n_values > 0) {
+            for (i = 0; i < n_values; i++) {
+                byte_pos += 4;      /* skip */
+                
+                len = *((int *)(byte_ptr + byte_pos));
+
+                byte_pos += 4;
+
+                byte_pos += 8;
+
+                if (i < vm && i < MAXVM) {
+                    value[i] = (byte_ptr + byte_pos);
+                }
+
+                byte_pos += ((len + 3) / 4) * 4;
+            }
+        }
+
+        if (!strcmp(name, "DiffusionDirectionality")) {
+            acr_insert_string(&group_list, ACR_Diffusion_directionality,
+                              value[0]);
+        }
+        else if (!strcmp(name, "B_value")) {
+            double tmp = atof(value[0]);
+            acr_insert_double(&group_list, ACR_Diffusion_b_value, 1, &tmp);
+        }
+        else if (!strcmp(name, "DiffusionGradientDirection")) {
+            double tmp[3];
+            if (vm == 3 && n_values >= vm) {
+                tmp[0] = atof(value[0]);
+                tmp[1] = atof(value[1]);
+                tmp[2] = atof(value[2]);
+                
+                acr_insert_double(&group_list, 
+                                  ACR_Diffusion_gradient_orientation,
+                                  3,
+                                  tmp);
+            }
+        }
+
+#if 0
+        if (G.Debug >= HI_LOGGING) {
+            printf("%s VM=%d VR=%s %d ", name, vm, vr, n_values);
+
+            if (n_values != 0) {
+                for (i = 0; i < vm && i < MAXVM; i++) {
+                    printf("%s ", value[i]);
+                }
+            }
+            printf("\n");
+        }
+#endif
+    }
+    return (group_list);
+}
+
 static Acr_Group 
 add_siemens_info(Acr_Group group_list)
 {
@@ -701,6 +827,11 @@
     char *str_ptr;
     int interpolation_flag;
 
+    element = acr_find_group_element(group_list, SPI_Protocol2);
+    if (element != NULL) {
+        group_list = parse_siemens_proto2(group_list, element);
+    }
+
     /* now fix the group list to provide essential info
      * via standard dicom elements
      * (this lets the rest of the code be more reusable)
@@ -725,6 +856,8 @@
 
     protocol = acr_find_group_element(group_list, SPI_Protocol);
     if (protocol != NULL) {
+        /* parse_siemens_junk(protocol); */
+
         if (G.Debug >= HI_LOGGING) {
             printf("Incorporating Siemens protocol structure...\n");
         }
@@ -901,6 +1034,56 @@
                 acr_insert_numeric(&group_list, ACR_Acquisition, enc_ix);
             }
         } /* end of diffusion scan handling */
+
+        /* There is another whole class of (probably newer) diffusion
+         * weighted images that use a very different arrangement and
+         * need better handling.
+         */
+        prot_find_string(protocol, "sDiffusion.ucDiffWeightedImage", str_buf);
+        if (!strcmp(str_buf, "0x1")) {
+            int enc_ix;
+
+            /* Flag this as a diffusion image.
+             */
+            acr_insert_string(&group_list, ACR_Acquisition_contrast, 
+                              "DIFFUSION");
+
+            prot_find_string(protocol, "sDiffusion.lDiffDirections", str_buf);
+            
+            acr_insert_numeric(&group_list, ACR_Acquisitions_in_series,
+                               atoi(str_buf) + 1);
+                
+            str_ptr = strchr(acr_find_string(group_list,
+                                             ACR_Sequence_name, ""), '#');
+            if (str_ptr == NULL) {
+                enc_ix = 0;
+            } 
+            else {
+                enc_ix = atoi(str_ptr + sizeof(char)) + 1;
+            }
+            acr_insert_numeric(&group_list, ACR_Acquisition, enc_ix);
+
+            /* BUG! TODO! FIXME!  In dcm2mnc.c the sequence name is
+             * used as one of the criteria for starting a new
+             * file. For a DTI sequence, we don't want this to
+             * happen. For now I replace the sequence name with a
+             * bogus value in order to keep the DTI scan from being
+             * split into pieces.  This isn't right. We should deal
+             * with this in a more graceful way.
+             */
+            acr_insert_string(&group_list, ACR_Sequence_name, "DTI");
+
+            /* Reset acquisition time to series time plus the series
+             * index.  Otherwise each slice may get its own
+             * time. This is also probably wrong and in the ideal
+             * world we should deal with this gracefully. But we have
+             * no good way of storing per-slice timing information
+             * right now.
+             */
+            acr_insert_numeric(&group_list, ACR_Acquisition_time,
+                               acr_find_double(group_list, ACR_Series_time, 0) 
+                               + enc_ix);
+        }
     }
     else {
         /* If no protocol dump was found we have to assume no
@@ -1128,6 +1311,7 @@
     char *str_ptr;
     int slice_count;
     int slice_index;
+    char str_buf[128];
 
     /* To use the Philips proprietary group, we have to figure out the
      * DICOM private creator ID in use.  The group ID is always 0x2001,
@@ -1158,8 +1342,88 @@
             }
         }
     }
-    if (creator_id.element_id > 0xff) {
+    else {
+        creator_id.element_id = 0;
+    }
+    if (creator_id.element_id > 0xff || creator_id.element_id < 0x01) {
         printf("WARNING: Can't find Philips private creator ID.\n");
+
+        /* OK, this may be an old Philips file with the SPI stuff in it.
+         */
+        str_ptr = acr_find_string(group_list, SPI_PMS_grp19_tag, "");
+        if (!strncmp(str_ptr, "PHILIPS MR", 10)) {
+            if (G.Debug >= HI_LOGGING) {
+                printf("Found Philips SPI information\n");
+            }
+
+            str_ptr = acr_find_string(group_list, ACR_Image_position_patient, "");
+            if (*str_ptr == '\0') {
+                double x, y, z;
+
+                x = acr_find_double(group_list, SPI_PMS_lr_position2, 0);
+                y = acr_find_double(group_list, SPI_PMS_ap_position2, 0);
+                z = acr_find_double(group_list, SPI_PMS_cc_position2, 0);
+                sprintf(str_buf, "%.15g\\%.15g\\%.15g", x, y, z);
+                printf("New position: %s\n", str_buf);
+                acr_insert_string(&group_list, ACR_Image_position_patient, 
+                                  str_buf);
+
+            }
+            str_ptr = acr_find_string(group_list, ACR_Number_of_slices, "");
+            if (*str_ptr == '\0') {
+                short n;
+                n = acr_find_short(group_list, SPI_PMS_slice_count, 0);
+                acr_insert_short(&group_list, ACR_Number_of_slices, n);
+            }
+            str_ptr = acr_find_string(group_list, 
+                                      ACR_Image_orientation_patient, "");
+            if (strchr(str_ptr, '\\') == NULL) {
+                int orientation = acr_find_int(group_list, 
+                                               SPI_PMS_slice_orientation,
+                                               1);
+                switch (orientation) {
+                case 1:         /* transverse, slice=Z */
+                    strcpy(str_buf, "1\\0\\0\\0\\1\\0");
+                    break;
+                case 2:         /* sagittal, slice=X */
+                    strcpy(str_buf, "0\\1\\0\\0\\0\\1");
+                    break;
+                case 3:         /* coronal, slice=Y */
+                    strcpy(str_buf, "1\\0\\0\\0\\0\\1");
+                    break;
+                }
+
+                acr_insert_string(&group_list, ACR_Image_orientation_patient,
+                                  str_buf);
+            }
+#if 0 /* not clear that this is needed */
+            str_ptr = acr_find_string(group_list, ACR_Pixel_size, "");
+            if (*str_ptr == '\0') {
+                pms_element = acr_find_group_element(group_list, 
+                                                     SPI_PMS_field_of_view);
+                if (pms_element != NULL) {
+                    double fov[2];
+                    acr_get_element_numeric_array(pms_element, 2, fov);
+                    if (fov[0] != 0.0 && fov[1] != 0.0) {
+                        double derived_spacing[2];
+                        int rows = acr_find_int(group_list, ACR_Rows, 1);
+                        int cols = acr_find_int(group_list, ACR_Columns, 1);
+
+                        derived_spacing[0] = fov[0] / cols;
+                        derived_spacing[1] = fov[1] / rows;
+
+                        sprintf(str_buf, "%.15g\\%.15g", 
+                                derived_spacing[0], derived_spacing[1]);
+
+                        printf("PHILIPS: New pixel size %s\n", str_buf);
+
+                        acr_insert_string(&group_list, ACR_Pixel_size, 
+                                          str_buf);
+                    }
+                }
+            }
+#endif
+        }
     }
     else {
         if (G.Debug >= HI_LOGGING) {