/* ---------------------------------------------------------------------- * code that deals with the CSA(1/2) header in Siemens mosaic images * * R Reynolds May 2011 * * see http://nipy.sourceforge.net/nibabel/dicom/siemens_csa.html * (which seems to be based on SPMs implementation) * * select notes from that page: * * - the field contents are always little-endian * - the CSA Image Header Info (0029 1010) field can be of two * types, CSA1 and CSA2 * * CSA1 format: * * 0029 1010 : S8, 8 byte field type * header: * uint32 n_tags : number of tags, in [1,128] else fail * uint32 unused : val == 77 * each tag: * S64 name : nul-terminated name string * int32 vm : ?? * S4 vr : possibly nul-term string * int32 syngodt : ?? * int32 nitems : number of items in CSA * int32 xx : maybe == 77 or 205 * each item: * int32 xx[4] : length (in bytes, maybe modified), ??, ??, ?? * SPM: if xx[0]-nitems < 0 or > remain bytes: * then terminate * uint8 value[item_len] : actual data, item_len = xx[0]-nitems * * CSA2 format: * * 0029 1010 : S8, 8 byte field type * header: * S4 hdr_id : val == "SV10" * uint8 unused1[4] : * uint32 n_tags : number of tags (in [1,128], else fail) * uint32 unused2 : val == 77 * each tag (same as for CSA1): * S64 name : nul-terminated name string * int32 vm : ?? * S4 vr : possibly nul-term string * int32 syngodt : ?? * int32 nitems : number of items in CSA * int32 xx : maybe == 77 or 205 * each item: * int32 xx[4] : length (in bytes), ??, ??, ?? * Note: interpretation differs from above * if xx[0] < 0 or > remain bytes, terminate * uint8 value[xx[0]] : actual data (length does not used nitems) * * ---------------------------------------------------------------------- */ typedef struct { unsigned int ntags; /* in [1,128], tags in CSA */ unsigned int unused; /* should be 77 */ } siemens_csa1_header; typedef struct { char hdr_id[4]; /* should be "SV10" */ unsigned char unused1[4]; /* in [1,128] */ unsigned int ntags; /* in [1,128], tags in CSA */ unsigned int unused2; /* should be 77 */ } siemens_csa2_header; typedef struct { char name[64]; /* nul-terminated string */ int vm; /* ?? */ char vr[4]; /* first 3 chars usable */ int syngodt; /* ?? */ int nitems; /* items under current tag */ int xx; /* 77 or 205 */ } siemens_csa_tag; typedef struct { int xx[4]; /* length, ??, ??, ?? */ unsigned char value[1]; /* data (len or len-nitems) */ } siemens_csa_item; static void swap_4(char * data, int nquads) { char * dp = data; char c, index; for( index = 0; index < nquads; index++ ) { c = dp[0]; dp[0] = dp[3]; dp[3] = c; c = dp[1]; dp[1] = dp[2]; dp[2] = c; dp += 4; } } static void swap_csa1_header(siemens_csa1_header * hdr) { swap_4((char *)&hdr->ntags, 1); swap_4((char *)&hdr->unused, 1); } static void swap_csa2_header(siemens_csa2_header * hdr) { swap_4((char *)&hdr->ntags, 1); swap_4((char *)&hdr->unused2, 1); } static void swap_csa_tag(siemens_csa_tag * hdr) { swap_4((char *)&hdr->vm, 1); swap_4((char *)&hdr->syngodt, 1); swap_4((char *)&hdr->nitems, 1); swap_4((char *)&hdr->xx, 1); } static void swap_csa_item(siemens_csa_item * hdr) { swap_4((char *)&hdr->xx, 4); } static int little_endian(void) { int val = 1; char * cp = (char *)&val; return *cp; } /* return 1 or 2, for CSA1 or CSA2 */ /* (CSA2 starts with SV10) */ static int get_csa_type(unsigned char * str) { if( !strncmp((char *)str, "SV10", 4) ) return 2; else return 1; } static int print_csa1_header(siemens_csa1_header * header, FILE * stream) { int size = sizeof(siemens_csa1_header); fprintf(stream, "-- CSA 1 header (size = %d)\n", size); fprintf(stream, " ntags : %u (in [1,128]?)\n", header->ntags); fprintf(stream, " unused : %u (77?)\n", header->unused); fputc('\n', stream); return size; } static int print_csa2_header(siemens_csa2_header * hdr, FILE * stream) { int size = sizeof(siemens_csa2_header); fprintf(stream, "-- CSA 2 header (size = %d)\n", size); fprintf(stream, " hdr_id : %4.4s (SV10?)\n", hdr->hdr_id); fprintf(stream, " unused1 : %u %u %u %u\n", hdr->unused1[0], hdr->unused1[1], hdr->unused1[2], hdr->unused1[3] ); fprintf(stream, " ntags : %u (in [1,128]?)\n", hdr->ntags); fprintf(stream, " unused2 : %u (77?)\n", hdr->unused2); fputc('\n', stream); return size; } static int print_csa_tag(siemens_csa_tag * tag, FILE * stream) { int size = sizeof(siemens_csa_tag); fprintf(stream, "-- CSA tag (size = %d)\n", size); fprintf(stream," tag name : %-64.64s\n", tag->name); fprintf(stream," vm : %d\n", tag->vm); fprintf(stream," vr : %-4.4s\n", tag->vr); fprintf(stream," syngodt : %d\n", tag->syngodt); fprintf(stream," nitems : %d\n", tag->nitems); fprintf(stream," xx : %d (77 or 205?)\n", tag->xx); fputc('\n', stream); return size; } /* print 'mesg' and then 'data' of length 'len' to 'fp' * if 'text': print as characters */ static int print_raw_data(char * mesg, unsigned char * data, int len, FILE * fp, int text) { int mlen = 0, posn = 0, indent = 0, nbytes; if( mesg ) { mlen = strlen(mesg); fputs(mesg, fp); } while( posn < len ) { if( posn % 32 == 0 ) { if( mesg ) { /* then indent, unless this is the first time */ if( indent > 0 ) fprintf(fp, "\n%*s", indent, ""); else indent = mlen; } } if( text && isprint(data[posn]) ) fputc(data[posn], fp); else if( text ) fprintf(fp, " 0x%02x", data[posn]); else fprintf(fp, " %02x", data[posn]); posn++; } fputc('\n', fp); fputc('\n', fp); return 0; } /* return xx[0] rounded up to multiple of 4 */ static int csa_item_value_size_aligned(siemens_csa_item * citem, int ilenoff) { int nbytes = citem->xx[0], nbytes_plus; int total_size; if( nbytes % 4 ) return nbytes + 4 - (nbytes % 4); else return nbytes; } static int csa_item_size(siemens_csa_item * citem, int ilenoff) { int size = 4*sizeof(int); /* base offset */ int nbytes = citem->xx[0], nbytes_plus; int total_size; total_size = size - ilenoff + csa_item_value_size_aligned(citem, ilenoff); return total_size; } /* print presumably nul-terminated string */ static int print_csa_item_str(siemens_csa_item * citem, int ilenoff, FILE *fp) { int nbytes = citem->xx[0]; if( nbytes > 0 ) fprintf(fp, " %.*s", nbytes, citem->value); return nbytes ; } /* - subtract ilenoff from length * - if data, print data (as text?) * - if remain is too small, terminate */ static int print_csa_item(siemens_csa_item * citem, int ilenoff, int data, int remain, FILE * stream) { int size = 4*sizeof(int); /* base offset */ int nbytes = citem->xx[0], nbytes_plus; int total_size; /* round nbytes up to next multiple of 4 */ nbytes_plus = csa_item_value_size_aligned(citem, ilenoff); total_size = csa_item_size(citem, ilenoff); fprintf(stream, "-- CSA item (size = %d, total = %d, remain = %d)\n", size, total_size, remain); fprintf(stream, " xx : %d (offset %d)\n", nbytes, ilenoff); fputc('\n', stream); if( total_size > remain ) { if(g_MDH_verb > 1) fprintf(stderr,"** PCSAI: insufficient space\n\n"); return total_size; } if( data && nbytes > 0 ) print_raw_data(" item data : ", citem->value, nbytes, stderr, 1); return total_size; } static int insert_slice_time(float st) { siemens_slice_times_t * ST = & g_siemens_slice_times; /* global struct */ ST->nused++; if( ST->nused > ST->nalloc ) { ST->nalloc = ST->nused + 4; /* alloc a few extra */ ST->times = (float *)realloc(ST->times, ST->nalloc*sizeof(float )); if( ! ST->times ) { fprintf(stderr,"** MIST: failed malloc of %d floats\n", ST->nalloc); ST->nused = 0; ST->nalloc = 0; return 1; } } ST->times[ST->nused-1] = st; return 0; } /* loop through all CSA data fields * * if name, only print that field contents * if data, print item data (1=raw, 2=text) * verb is verbose level * */ static int process_csa_data(unsigned char * str, int len, int verb, int data, char * name) { siemens_csa1_header * c1ptr = (siemens_csa1_header *)str; siemens_csa2_header * c2ptr = (siemens_csa2_header *)str; siemens_csa_tag * ctag = NULL; siemens_csa_item * citem = NULL; unsigned char * ucp = (unsigned char *)str; /* main pointer */ unsigned int itag, ntags; char * endp; float stime = 0.0; int tagsize = sizeof(siemens_csa_tag); int posn, ctype, remain, nitem, iitem, ilen, ilenoff=-1; int little = little_endian(); if( !str ) { if(verb > 1) fprintf(stderr,"** PACSAD: bad data pointer\n"); return 1; } if( verb > 2 ) { if( little ) fprintf(stderr,"-- little endian, so no swapping\n"); else fprintf(stderr,"-- big endian, will swap structures\n"); } ctype = get_csa_type(ucp); posn = 0; /* current index into ucp/str */ if( ctype == 1 ) { if( ! little ) swap_csa1_header(c1ptr); if( verb > 1 ) print_csa1_header(c1ptr, stderr); posn += sizeof(siemens_csa1_header); ntags = c1ptr->ntags; } else { if( ! little ) swap_csa2_header(c2ptr); if( verb > 1 ) print_csa2_header(c2ptr, stderr); posn += sizeof(siemens_csa2_header); ntags = c2ptr->ntags; } if( ntags > 128 ) { if(verb > 1) fprintf(stderr,"** PACSAD: ntags (%d) > 128\n", ntags); return 1; } for( itag = 0; itag < ntags; itag++ ) { remain = len - posn; if( remain < tagsize ) { if(verb>1) fprintf(stderr,"** PACSAD: no room for next tag\n"); return 1; } ctag = (siemens_csa_tag *)(ucp+posn); if( ! little ) swap_csa_tag(ctag); if( data || verb > 2 ) print_csa_tag(ctag, stderr); /* check for our name of interest */ if ( name && !strcmp(ctag->name, name) ) { if( verb > 1 ) print_csa_tag(ctag, stderr); g_siemens_slice_times.nused = 0; /* clear slice times */ } posn += sizeof(siemens_csa_tag); /* if CSA1 and FIRST tag, note the item length offset */ if( ilenoff < 0 ) { if( ctype == 1 ) { ilenoff = ctag->nitems; if(verb > 2) fprintf(stderr,"-- CSA1 item offset = %d\n",ilenoff); if( ilenoff < 0 ) ilenoff = 0; } else ilenoff = 0; } for( iitem = 0; iitem < ctag->nitems; iitem++ ) { citem = (siemens_csa_item *)(ucp+posn); if( ! little ) swap_csa_item(citem); if ( name && !strcmp(ctag->name, name) ) { if(verb==3) print_csa_item(citem, ilenoff, data, len-posn, stderr); if( verb > 1 ) print_csa_item_str(citem, ilenoff, stderr); /* our main purpose, actually insert each slice time */ stime = strtod((char *)citem->value, &endp); if( endp && endp > (char *)citem->value ) insert_slice_time(stime); else if (verb > 2) fprintf(stderr," "); } if( data || verb > 3 ) print_csa_item(citem, ilenoff, data, len-posn, stderr); posn += csa_item_size(citem, ilenoff); if( posn > len ) break; } /* if we have processed our tag of interest, we're done */ if ( name && !strcmp(ctag->name, name) ) { if( verb > 1 ) fputc('\n', stderr); if( verb < 5 ) return 0; /* else print everything */ } } return 0; } /* quick strstr function, but limited by len (not NUL char) */ static char * findstr(char * instr, char * sstr, int len) { char * cptr = instr; int posn = 0, slen; if( !instr || !sstr || len <= 0 ) return NULL; slen = strlen(sstr); len -= slen; while( posn <= len ){ if( *cptr != *sstr ) { cptr++; posn++; continue; } if( ! strncmp(cptr, sstr, slen) ) return cptr; else { cptr++; posn++; } } return NULL; } /* - search field '0x0029 1010' for string MosaicRefAcqTimes * - from there, find text formatted floats * - stop at AutoInlineImageFilterEnabled, if found */ static int check_for_mosaic_slice_times(PRV_ELEMENT_ITEM * elementItem) { unsigned el_gr = DCM_TAG_GROUP(elementItem->element.tag); unsigned el_el = DCM_TAG_ELEMENT(elementItem->element.tag); int el_len = elementItem->element.length; char start_txt[] = "MosaicRefAcqTimes"; char end_txt[] = "AutoInlineImageFilterEnabled"; siemens_slice_times_t * ST = & g_siemens_slice_times; /* global struct */ char * instr, * mstr; /* input string and Mosaic string addr */ char * s2; /* second search string, posn of AutoInline... */ char * pstr; /* position pointer, for reading times */ int rem, rem2 = 0; /* remainder counts */ int off, c, rv, diff; /* offset and counter vars */ float stime; /* any read slice time */ /* if this is not the correct element, nothing to do */ if( el_gr != 0x0029 || el_el != 0x1010 ) return 0; /* we are in the correct field, start by clearing old results */ ST->nused = 0; /* input string is field text, mstr is resulting MosaicRef text pointer */ instr = (char *)elementItem->element.d.ot; mstr = findstr(instr, start_txt, el_len); if( ! mstr ) { if( g_MDH_verb > 2 ) fprintf(stderr, "-- CFMST, no Mosaic string\n"); return 0; } off = mstr - instr; /* offset of Mosaic string into field */ rem = el_len - off; /* remaining length of field */ /* secondary remainder to be until any AutoInline... string */ s2 = findstr(mstr, end_txt, rem); if( s2 ) rem2 = s2 - mstr; else rem2 = 0; if( g_MDH_verb > 1 ) fprintf(stderr, "== found %s in 0x0029 1010\n" " off = %d of %d bytes (rem %d, %d)\n", start_txt, off, el_len, rem, rem2); if( s2 ) rem = rem2; /* after verbose, update remaining length */ if( rem <= 0 ) return 0; process_csa_data((unsigned char *)instr, el_len, g_MDH_verb, g_MDH_verb>3, start_txt); /* in really verbose mode, print out raw text */ if( g_MDH_verb > 3 ) { unsigned char * ucp = (unsigned char *)mstr; fprintf(stderr, "-- remaining spaced hex or digit or '.' :\n"); for(c=0; c 1 ) fprintf(stderr,"\n++ found %d slice times\n", ST->nused); return 0; }