~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina4/src/file/dicom/dicom_io.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /**@(#)Functions to read/write DICOM image files
  2        These functions are a modified version of the 
  3        ACR-NEMA ones
  4 */
  5 
  6 #include <stdio.h>
  7 #include <string.h>
  8 #include <tina/file.h>
  9 #include <tina/sys.h>
 10 #include <tina/sys_types.h>
 11 #include <tina/sysfuncs.h>
 12 #include <tina/math.h>
 13 #include <tina/mathfuncs.h>
 14 #include <tina/temporal.h>
 15 #include <tina/dicom_dic.h>
 16 
 17 #define VOXELS    450
 18 #define DYNSTIME  451
 19 #define TE_DATA   453
 20 #define PAT_DATA  454
 21 
 22 #define HEADERMAXCOUNT 100
 23 #define PART10    1
 24 #define NONPART10 0
 25 #define RUBBISH   -1
 26 
 27 #define I_LITTLE_ENDIAN "1.2.840.10008.1.2"
 28 #define E_LITTLE_ENDIAN "1.2.840.10008.1.2.1"
 29 #define E_BIG_ENDIAN    "1.2.840.10008.1.2.2"
 30 
 31 
 32 extern Bool fclose_2(FILE * stream, const char *pathname);
 33 extern FILE *fopen_2(const char *pathname, const char *mode);
 34 extern Bool fread_imrect_data(const Imrect * imrect, FILE * stream,
 35                                                                                                                         const char *pathname);
 36 extern Bool fwrite_imrect_data(const Imrect * imrect, FILE * stream,
 37                                                                                                                          const char *pathname);
 38 extern void dynstimes_free(float *times);
 39 extern void im_endian_inplace(Imrect * imrect);
 40 
 41 static int dicom_format = RUBBISH;
 42 
 43 
 44 static void set_dicom_format(int format)
 45 {
 46         dicom_format = format;
 47 }
 48 
 49 
 50 static int dblock_vr_conv(unsigned int *dblock, FILE * fp)
 51 {
 52         size_t readerr;
 53         char val, junk[32];
 54         unsigned int   dtemp;
 55         unsigned char *pvr;
 56 
 57         if (dicom_format == RUBBISH)
 58                 return (-1);
 59 
 60         pvr = (char *) dblock;
 61 
 62         /* 
 63          * VR is a nested sequence (SQ)
 64          */
 65 
 66         if ((pvr[0] == 'S' && pvr[1] == 'Q') || (pvr[3] == 'S' && pvr[2] == 'Q'))
 67         {
 68                 readerr = fread(dblock, sizeof(int), 1, fp);
 69                 word_swap((char *) dblock);
 70                 *dblock = 0;
 71                 return (1);
 72         }
 73 
 74         if ((pvr[0] == 'O' && pvr[1] == 'B') || (pvr[0] == 'O' && pvr[1] == 'W'))
 75         {
 76                 readerr = fread(dblock, sizeof(int), 1, fp);
 77                 word_swap((char *) dblock);
 78                 return (1);
 79         }
 80 
 81         if ((pvr[3] == 'O' && pvr[2] == 'B') || (pvr[3] == 'O' && pvr[2] == 'W'))
 82         {
 83                 readerr = fread(dblock, sizeof(int), 1, fp);
 84                 word_swap((char *) dblock);
 85                 return (1);
 86         }
 87 
 88         if (*dblock == 0xffffffff)
 89         {
 90                 *dblock = 0;
 91                 return (1);
 92         }
 93 
 94         if ((pvr[0] == 'A' && pvr[1] == 'E') ||
 95                         (pvr[0] == 'A' && pvr[1] == 'S') ||
 96                         (pvr[0] == 'A' && pvr[1] == 'T') ||
 97                         (pvr[0] == 'C' && pvr[1] == 'S') ||
 98                         (pvr[0] == 'D' && pvr[1] == 'A') ||
 99                         (pvr[0] == 'D' && pvr[1] == 'S') ||
100                         (pvr[0] == 'D' && pvr[1] == 'T') ||
101                         (pvr[0] == 'F' && pvr[1] == 'L') ||
102                         (pvr[0] == 'F' && pvr[1] == 'D') ||
103                         (pvr[0] == 'I' && pvr[1] == 'S') ||
104                         (pvr[0] == 'L' && pvr[1] == 'O') ||
105                         (pvr[0] == 'L' && pvr[1] == 'T') ||
106                         (pvr[0] == 'P' && pvr[1] == 'N') ||
107                         (pvr[0] == 'S' && pvr[1] == 'H') ||
108                         (pvr[0] == 'S' && pvr[1] == 'L') ||
109                         (pvr[0] == 'S' && pvr[1] == 'S') ||
110                         (pvr[0] == 'S' && pvr[1] == 'T') ||
111                         (pvr[0] == 'T' && pvr[1] == 'M') ||
112                         (pvr[0] == 'U' && pvr[1] == 'I') ||
113                         (pvr[0] == 'U' && pvr[1] == 'L') || (pvr[0] == 'U' && pvr[1] == 'S'))
114         {
115                 /*
116                 (*dblock) &= 0xffff0000;
117                 (*dblock) >>= 16;
118                 */
119 
120 #ifdef LITTLE_ENDIAN_ARCHITECTURE
121                 dtemp = (unsigned int)(256*pvr[3] + pvr[2]);
122 #else
123                 dtemp = (unsigned int)(256*pvr[2] + pvr[3]);
124 #endif
125 
126                 (*dblock) = dtemp;
127                 return (1);
128         }
129 
130         if ((pvr[3] == 'A' && pvr[2] == 'E') ||
131                         (pvr[3] == 'A' && pvr[2] == 'S') ||
132                         (pvr[3] == 'A' && pvr[2] == 'T') ||
133                         (pvr[3] == 'C' && pvr[2] == 'S') ||
134                         (pvr[3] == 'D' && pvr[2] == 'A') ||
135                         (pvr[3] == 'D' && pvr[2] == 'S') ||
136                         (pvr[3] == 'D' && pvr[2] == 'T') ||
137                         (pvr[3] == 'F' && pvr[2] == 'L') ||
138                         (pvr[3] == 'F' && pvr[2] == 'D') ||
139                         (pvr[3] == 'I' && pvr[2] == 'S') ||
140                         (pvr[3] == 'L' && pvr[2] == 'O') ||
141                         (pvr[3] == 'L' && pvr[2] == 'T') ||
142                         (pvr[3] == 'P' && pvr[2] == 'N') ||
143                         (pvr[3] == 'S' && pvr[2] == 'H') ||
144                         (pvr[3] == 'S' && pvr[2] == 'L') ||
145                         (pvr[3] == 'S' && pvr[2] == 'S') ||
146                         (pvr[3] == 'S' && pvr[2] == 'T') ||
147                         (pvr[3] == 'T' && pvr[2] == 'M') ||
148                         (pvr[3] == 'U' && pvr[2] == 'I') ||
149                         (pvr[3] == 'U' && pvr[2] == 'L') || (pvr[3] == 'U' && pvr[2] == 'S'))
150         {
151                 /*
152                 (*dblock) &= 0x0000ffff;
153                 */
154 
155 #ifdef LITTLE_ENDIAN_ARCHITECTURE
156                 dtemp = (unsigned int)(256*pvr[1] + pvr[0]);
157 #else
158                 dtemp = (unsigned int)(256*pvr[0] + pvr[1]);
159 #endif
160 
161                 (*dblock) =  dtemp;
162                 return (1);
163         }
164 
165         /*
166                  must be implicit VR
167   */
168         return (-1);
169 }
170 
171 /*
172   finds endian UID in part10
173   leaves file ptr after last DCM_GROUPFILEMETA tag
174 */
175 static Bool dicom_part10_endian_swap(FILE * fp)
176 {
177         Bool success = false;
178         char endian[32];
179         unsigned char junk;
180         unsigned short group, element;
181         unsigned int i, dblock;
182         long int where = ftell(fp);
183 
184 /* part10 meta dicom are implicit little endian by default */
185 
186 #ifdef LITTLE_ENDIAN_ARCHITECTURE
187         Bool little = false;
188         Bool big = true;
189 
190         set_swapping_ts(0);
191 #else
192         Bool little = true;
193         Bool big = false;
194 
195         set_swapping_ts(1);
196 #endif
197 
198         fread(&group, sizeof(short), 1, fp);
199         short_swap((char *) &group);
200 
201         while (fp && group <= DCM_GROUPFILEMETA)
202         {
203                 fread(&element, sizeof(short), 1, fp);
204                 short_swap((char *) &element);
205                 fread(&dblock, sizeof(int), 1, fp);
206                 word_swap((char *) &dblock);
207                 dblock_vr_conv(&dblock, fp);
208 
209                 switch (DCM_MAKETAG(group, element))
210                 {
211                         case DCM_METATRANSFERSYNTAX:
212                                 fread(endian, sizeof(char), dblock, fp);
213                                 endian[dblock] = '\0';
214                                 success = true;
215                                 break;
216 
217                         case DCM_DLMITEM:
218                         case DCM_DLMITEMDELIMITATIONITEM:
219                         case DCM_DLMSEQUENCEDELIMITATIONITEM:
220                                 break;
221 
222                         default:
223                                 for (i = 0; i < dblock; i++)
224                                         fread(&junk, sizeof(char), 1, fp);
225                                 break;
226                 }
227 
228                 where = ftell(fp);
229                 fread(&group, sizeof(short), 1, fp);
230                 short_swap((char *) &group);
231         }
232 
233 /* put fptr back before last group read */
234 
235         fseek(fp, where, SEEK_SET);
236 
237         if (success && !strcmp(endian, E_BIG_ENDIAN))
238                 return (big);
239 
240         return (little);
241 }
242 
243 
244 /*
245   part 10 dicom has a preamble of 128 bytes followed by
246   "DICM", although preamble may be missing. 
247   sets endian.
248 */
249 static int dicom_preread_tests(FILE * fp)
250 {
251         size_t readerr;
252         short group, element, temp;
253         char td, ti, tc, tm;
254         int i, dblock, repeat, count;
255 
256         if (fp == NULL)
257                 return (RUBBISH);
258 
259 /* part10 DICM header */
260         fseek(fp, 128, SEEK_SET);
261         readerr = fread(&td, sizeof(char), 1, fp);
262         readerr = fread(&ti, sizeof(char), 1, fp);
263         readerr = fread(&tc, sizeof(char), 1, fp);
264         readerr = fread(&tm, sizeof(char), 1, fp);
265         if (td == 'D' && ti == 'I' && tc == 'C' && tm == 'M')
266         {
267                 set_dicom_format(PART10);
268 
269                 if (dicom_part10_endian_swap(fp))
270                         set_swapping_ts(1);
271                 else
272                         set_swapping_ts(0);
273 
274                 return (PART10);
275         }
276 
277 /* no pre-amble part10 DICM header */
278         rewind(fp);
279         readerr = fread(&td, sizeof(char), 1, fp);
280         readerr = fread(&ti, sizeof(char), 1, fp);
281         readerr = fread(&tc, sizeof(char), 1, fp);
282         readerr = fread(&tm, sizeof(char), 1, fp);
283         if (td == 'D' && ti == 'I' && tc == 'C' && tm == 'M')
284         {
285                 set_dicom_format(PART10);
286 
287                 if (dicom_part10_endian_swap(fp))
288                         set_swapping_ts(1);
289                 else
290                         set_swapping_ts(0);
291 
292                 return (PART10);
293         }
294 
295         repeat = 0;
296         set_swapping_ts(0);
297 
298         while (repeat < 2)
299         {
300                 rewind(fp);
301                 readerr = fread(&group, sizeof(short), 1, fp);
302                 short_swap((char *) &group);
303                 readerr = fread(&element, sizeof(short), 1, fp);
304                 short_swap((char *) &element);
305                 readerr = fread(&dblock, sizeof(int), 1, fp);
306                 word_swap((char *) &dblock);
307 
308                 count = 0;
309                 while (group < DCM_GROUPIDENTIFYING && count < HEADERMAXCOUNT)
310                 {
311                         for (i = 0; i < dblock; i++)
312                                 readerr = fread(&temp, sizeof(char), 1, fp);
313 
314                         readerr = fread(&group, sizeof(short), 1, fp);
315                         short_swap((char *) &group);
316                         readerr = fread(&element, sizeof(short), 1, fp);
317                         short_swap((char *) &element);
318                         readerr = fread(&dblock, sizeof(int), 1, fp);
319                         word_swap((char *) &dblock);
320                         count++;
321                 }
322 
323 /* standard non-part10 header start */
324                 if (DCM_MAKETAG(group, element) == DCM_IDGROUPLENGTH && dblock == 4)
325                 {
326                         rewind(fp);
327                         set_dicom_format(NONPART10);
328                         return (NONPART10);
329                 }
330 
331 /* inferior although common non-part10 header starts */
332          if ((DCM_MAKETAG(group, element) == DCM_IDLENGTHTOEND) ||
333              (DCM_MAKETAG(group, element) == DCM_IDIMAGETYPE) ||
334              (DCM_MAKETAG(group, element) == DCM_IDSPECIFICCHARACTER))
335                 {
336                         rewind(fp);
337                         set_dicom_format(NONPART10);
338                         return (NONPART10);
339                 }
340 
341                 set_swapping_ts(1);
342                 repeat++;
343         }
344 
345         return (RUBBISH);                                                       /* failed */
346 }
347 
348 
349 
350 Imrect *im_dicom_conv(Imrect * im)
351 {
352         Imrect *im2;
353         Imregion roi;
354         unsigned char *row1;
355         unsigned char *row2;
356         int i, j, k;
357         int lx, ux, ly, uy;
358 
359         if (im == NULL)
360                 return (NULL);
361 
362         roi = *(im->region);
363         roi.ux = roi.ux * 4.0 / 3;
364         im2 = im_alloc(im->height, roi.ux, &roi, short_v);
365         lx = roi.lx;
366         ux = roi.ux;
367         ly = roi.ly;
368         uy = roi.uy;
369 
370         for (i = ly; i < uy; ++i)
371         {
372                 row1 = IM_ROW(im, i);
373                 row2 = IM_ROW(im2, i);
374                 for (k = 2 * lx, j = 2 * lx; k < 2 * ux; j += 3)
375                 {
376                         row2[k++] = row1[j];
377                         row2[k++] = row1[j + 1] % 16;
378                         row2[k++] = row1[j + 1] / 16 + (row1[j + 2] % 16) * 16;
379                         row2[k++] = row1[j + 2] / 16;
380                 }
381         }
382 
383         return (im2);
384 }
385 
386 
387 Bool dicom_hdr_TE_extract(FILE * fp, Imrect * im)
388 {
389         Bool success = false;
390         float *TE;
391         int type;
392         int i, j, k;
393         char *buffer = NULL;
394         char temp[64];
395         unsigned char junk;
396         unsigned short group, element;
397         unsigned int dblock;
398 
399         if ((type = dicom_preread_tests(fp)) == RUBBISH)
400         {
401                 format("error: file failed pre checks\n");
402                 return (false);
403         }
404 
405         TE = (float *) ralloc(sizeof(float));
406         fread(&group, sizeof(short), 1, fp);
407         short_swap((char *) &group);
408 
409         while (fp && !success && group != DCM_GROUPPIXEL)
410         {
411                 fread(&element, sizeof(short), 1, fp);
412                 short_swap((char *) &element);
413                 fread(&dblock, sizeof(int), 1, fp);
414                 word_swap((char *) &dblock);
415                 dblock_vr_conv(&dblock, fp);
416 
417                 switch (DCM_MAKETAG(group, element))
418                 {
419                         case DCM_ACQECHOTIME:
420                                 fread(temp, sizeof(char), dblock, fp);
421                                 temp[dblock] = '\0';
422                                 sscanf(temp, "%f", TE);
423                                 *TE = *TE / 100.0;
424                                 success = true;
425                                 break;
426 
427                         case DCM_DLMITEM:
428                         case DCM_DLMITEMDELIMITATIONITEM:
429                         case DCM_DLMSEQUENCEDELIMITATIONITEM:
430                                 break;
431 
432                         default:
433                                 for (i = 0; i < dblock; i++)
434                                         fread(&junk, sizeof(char), 1, fp);
435                                 break;
436                 }
437                 fread(&group, sizeof(short), 1, fp);
438                 short_swap((char *) &group);
439         }
440 
441         if (success)
442         {
443                 im->props = proplist_rm(im->props, TE_DATA);
444                 im->props =
445                                 proplist_addifnp(im->props, (void *) TE, TE_DATA, rfree, false);
446         }
447 
448         return (success);
449 }
450 
451 
452 Bool dicom_hdr_dynstimes_extract(FILE * fp, Imrect * im)
453 {
454         Bool success = false;
455         float *times;
456         int type;
457         int i, j, k, stime;
458         int dynstudy;
459         char *buffer = NULL;
460         char temp[64], v;
461         unsigned char junk;
462         unsigned short group, element;
463         unsigned int dblock;
464 
465         if ((type = dicom_preread_tests(fp)) == RUBBISH)
466         {
467                 format("error: file failed pre checks\n");
468                 return (false);
469         }
470 
471         fread(&group, sizeof(short), 1, fp);
472         short_swap((char *) &group);
473 
474         while (fp && !success && group != DCM_GROUPPIXEL)
475         {
476                 fread(&element, sizeof(short), 1, fp);
477                 short_swap((char *) &element);
478                 fread(&dblock, sizeof(int), 1, fp);
479                 word_swap((char *) &dblock);
480                 dblock_vr_conv(&dblock, fp);
481 
482                 switch (DCM_MAKETAG(group, element))
483                 {
484                         case DCM_RELNUMBERTEMPORALPOSITIONS:
485                                 fread(temp, sizeof(char), dblock, fp);
486                                 temp[dblock] = '\0';
487                                 sscanf(temp, "%d", &stime);
488                                 break;
489 
490                         case DCM_MAKETAG(0x0019, 0x1021):
491                                 buffer = (char *) ralloc(dblock * sizeof(char));
492                                 times = fvector_alloc(0, stime);
493                                 fread(buffer, sizeof(char), dblock, fp);
494 
495                                 j = k = 0;
496                                 for (i = 0; i < stime; i++)
497                                 {
498                                         while (buffer[j] != '\\')
499                                                 temp[k++] = buffer[j++];
500                                         temp[k] = '\0';
501                                         sscanf(temp, "%f", &(times[i]));
502                                         k = 0;
503                                         j++;
504                                 }
505                                 rfree(buffer);
506                                 success = true;
507                                 break;
508 
509                         case DCM_DLMITEM:
510                         case DCM_DLMITEMDELIMITATIONITEM:
511                         case DCM_DLMSEQUENCEDELIMITATIONITEM:
512                                 break;
513 
514 
515                         default:
516                                 for (i = 0; i < dblock; i++)
517                                         fread(&junk, sizeof(char), 1, fp);
518                                 break;
519                 }
520                 fread(&group, sizeof(short), 1, fp);
521                 short_swap((char *) &group);
522         }
523 
524         if (success)
525         {
526                 im->props = proplist_rm(im->props, DYNSTIME);
527                 im->props =
528                                 proplist_addifnp(im->props, (void *) times, DYNSTIME,
529                                                                                                  dynstimes_free, false);
530         }
531 
532         return (success);
533 }
534 
535 
536 Bool dicom_hdr_heartrate_extract(FILE * fp, Imrect * im)
537 {
538         Bool success = false;
539         float *times, btime, phases, heartrate;
540         float beattime;
541         int type;
542         int i, j, k;
543         int psuccess = 0;
544         char *buffer = NULL;
545         char temp[64], v;
546         unsigned char junk;
547         unsigned short group, element;
548         unsigned int dblock;
549 
550         if ((type = dicom_preread_tests(fp)) == RUBBISH)
551         {
552                 format("error: file failed pre checks\n");
553                 return (false);
554         }
555 
556         fread(&group, sizeof(short), 1, fp);
557         short_swap((char *) &group);
558 
559         while (fp && !success && group != DCM_GROUPPIXEL)
560         {
561                 fread(&element, sizeof(short), 1, fp);
562                 short_swap((char *) &element);
563                 fread(&dblock, sizeof(int), 1, fp);
564                 word_swap((char *) &dblock);
565                 dblock_vr_conv(&dblock, fp);
566 
567                 switch (DCM_MAKETAG(group, element))
568                 {
569                         case DCM_MAKETAG(0x0019, 0x1022):
570                                 fread(temp, sizeof(char), dblock, fp);
571                                 temp[dblock] = '\0';
572                                 sscanf(temp, "%f", &btime);
573                                 if (psuccess == 2)
574                                         success = true;
575                                 else
576                                         psuccess++;
577                                 break;
578 
579                         case DCM_MAKETAG(0x0019, 0x1069):
580                                 fread(temp, sizeof(char), dblock, fp);
581                                 temp[dblock] = '\0';
582                                 sscanf(temp, "%f", &phases);
583                                 if (psuccess == 2)
584                                         success = true;
585                                 else
586                                         psuccess++;
587                                 break;
588 
589                         case DCM_MAKETAG(0x0019, 0x106A):
590                                 fread(temp, sizeof(char), dblock, fp);
591                                 temp[dblock] = '\0';
592                                 sscanf(temp, "%f", &heartrate);
593                                 if (psuccess == 2)
594                                         success = true;
595                                 else
596                                         psuccess++;
597                                 break;
598 
599                         case DCM_DLMITEM:
600                         case DCM_DLMITEMDELIMITATIONITEM:
601                         case DCM_DLMSEQUENCEDELIMITATIONITEM:
602                                 break;
603 
604 
605                         default:
606                                 for (i = 0; i < dblock; i++)
607                                         fread(&junk, sizeof(char), 1, fp);
608                                 break;
609                 }
610                 fread(&group, sizeof(short), 1, fp);
611                 short_swap((char *) &group);
612         }
613 
614         if (phases == 1)
615                 success = false;
616         if (success)
617         {
618                 times = fvector_alloc(0, phases);
619                 beattime = 1.0 / (heartrate / 60);
620                 for (i = 0; i < phases; i++)
621                         times[i] = btime + (float) i *(beattime / (phases - 1));
622                 im->props = proplist_rm(im->props, DYNSTIME);
623                 im->props =
624                                 proplist_addifnp(im->props, (void *) times, DYNSTIME, rfree,
625                                                                                                  false);
626         }
627 
628         return (success);
629 }
630 
631 
632 Bool dicom_hdr_patientdetails_extract(FILE * fp, Imrect * im)
633 {
634         Bool psuccess = false;
635         int i, type;
636         char *details = NULL;
637         unsigned char junk;
638         unsigned short group, element;
639         unsigned int dblock;
640 
641         if ((type = dicom_preread_tests(fp)) == RUBBISH)
642         {
643                 format("error: file failed pre checks\n");
644                 return (false);
645         }
646 
647         fread(&group, sizeof(short), 1, fp);
648         short_swap((char *) &group);
649 
650         details = (char *) cvector_alloc(0, 256);
651         while (fp && !psuccess && group != DCM_GROUPPIXEL)
652         {
653                 fread(&element, sizeof(short), 1, fp);
654                 short_swap((char *) &element);
655                 fread(&dblock, sizeof(int), 1, fp);
656                 word_swap((char *) &dblock);
657                 dblock_vr_conv(&dblock, fp);
658 
659                 switch (DCM_MAKETAG(group, element))
660                 {
661                         case DCM_MAKETAG(0x0008, 0x0020):
662                                 fread(details, sizeof(char), dblock, fp);
663                                 details[dblock] = '\0';
664                                 psuccess = true;
665                                 break;
666 
667                         case DCM_DLMITEM:
668                         case DCM_DLMITEMDELIMITATIONITEM:
669                         case DCM_DLMSEQUENCEDELIMITATIONITEM:
670                                 break;
671 
672 
673                         default:
674                                 for (i = 0; i < dblock; i++)
675                                         fread(&junk, sizeof(char), 1, fp);
676                                 break;
677                 }
678                 fread(&group, sizeof(short), 1, fp);
679                 short_swap((char *) &group);
680         }
681 
682         if (psuccess != true)
683                 cvector_free(details, 0);
684         else
685         {
686                 im->props = proplist_rm(im->props, PAT_DATA);
687                 im->props = proplist_addifnp(im->props, (void *) details,
688                                                                                                                                  PAT_DATA, rfree, false);
689         }
690 
691         return (psuccess);
692 }
693 
694 
695 Bool dicom_hdr_voxelscale_extract(FILE * fp, Imrect * im)
696 {
697         Vec3 *iscale;
698         float xsize, ysize, zsize;
699         int type;
700         int psuccess = 0;
701         int i;
702         char dimension[64];
703         unsigned char junk;
704         unsigned short group, element;
705         unsigned int dblock;
706 
707         xsize = ysize = zsize = 1.0;
708         if ((type = dicom_preread_tests(fp)) == RUBBISH)
709         {
710                 format("error: file failed pre checks\n");
711                 return (false);
712         }
713 
714         fread(&group, sizeof(short), 1, fp);
715         short_swap((char *) &group);
716 
717         while (fp && group != DCM_GROUPPIXEL && psuccess < 3)
718         {
719 
720                 fread(&element, sizeof(short), 1, fp);
721                 short_swap((char *) &element);
722                 fread(&dblock, sizeof(int), 1, fp);
723                 word_swap((char *) &dblock);
724                 dblock_vr_conv(&dblock, fp);
725 
726                 switch (DCM_MAKETAG(group, element))
727                 {
728                         case DCM_ACQSLICESPACING:
729                                 fread(dimension, sizeof(char), dblock, fp);
730                                 dimension[dblock] = '\0';
731                                 sscanf(dimension, "%f", &zsize);
732                                 psuccess++;
733                                 break;
734 
735                         case DCM_IMGPIXELSPACING:
736                                 fread(dimension, sizeof(char), dblock, fp);
737                                 dimension[dblock] = '\0';
738                                 sscanf(dimension, "%f\\%f", &xsize, &ysize);
739                                 psuccess++;
740                                 break;
741 
742                         case DCM_DLMITEM:
743                         case DCM_DLMITEMDELIMITATIONITEM:
744                         case DCM_DLMSEQUENCEDELIMITATIONITEM:
745                                 break;
746 
747 
748                         default:
749                                 for (i = 0; i < dblock; i++)
750                                         fread(&junk, sizeof(char), 1, fp);
751                                 break;
752                 }
753                 fread(&group, sizeof(short), 1, fp);
754                 short_swap((char *) &group);
755         }
756 
757         if (psuccess == 2)
758         {
759                 iscale = vec3_alloc();
760                 iscale->el[0] = xsize;
761                 iscale->el[1] = ysize;
762                 iscale->el[2] = zsize;
763 
764                 im->props = proplist_rm(im->props, VOXELS);
765                 im->props = proplist_addifnp(im->props, iscale,
766                                                                                                                                  VOXELS, vec3_free, false);
767         }
768 
769         return ((Bool) (psuccess == 3));
770 }
771 
772 static Imrect *dicom_read_multiformat_image(const char *pathname,
773                                                                                                                                                                                 FILE * fp_img)
774 {
775         int get_swapping_ts();
776         Imrect *imrect = NULL, *imrect2 = NULL;
777         Imregion imregion;
778         Vartype new_vtype;
779         float range;
780         float rinter, rslope;
781         float val, pix;
782         int i, j, k;
783         unsigned char junk;
784         char dimension[64];
785         unsigned short abits;
786         unsigned short sbits;
787         unsigned short sign;
788         unsigned short rows;
789         unsigned short cols;
790         unsigned short group;
791         unsigned short element;
792         unsigned int dblock;
793         float xsize, ysize, zsize;
794         Vec3 *iscale;
795 
796         rinter = 0.0;
797         rslope = 1.0;
798         xsize = ysize = zsize = 1.0;
799 
800         fread(&group, sizeof(short), 1, fp_img);
801         short_swap((char *) &group);
802         fread(&element, sizeof(short), 1, fp_img);
803         short_swap((char *) &element);
804 
805         while (fp_img && DCM_MAKETAG(group, element) != DCM_PXLPIXELDATA)
806         {
807                 fread(&dblock, sizeof(int), 1, fp_img);
808                 word_swap((char *) &dblock);
809                 dblock_vr_conv(&dblock, fp_img);
810 
811                 if (DCM_MAKETAG(group, element) == DCM_PXLPIXELDATA)
812                         break;
813 
814                 switch (DCM_MAKETAG(group, element))
815                 {
816                         case DCM_ACQSLICESPACING:
817                                 fread(dimension, sizeof(char), dblock, fp_img);
818                                 dimension[dblock] = '\0';
819                                 sscanf(dimension, "%f", &zsize);
820                                 break;
821 
822                         case DCM_IMGPIXELSPACING:
823                                 fread(dimension, sizeof(char), dblock, fp_img);
824                                 dimension[dblock] = '\0';
825                                 sscanf(dimension, "%f\\%f", &xsize, &ysize);
826                                 break;
827 
828                         case DCM_IMGRESCALEINTERCEPT:
829                                 fread(dimension, sizeof(char), dblock, fp_img);
830                                 dimension[dblock] = '\0';
831                                 sscanf(dimension, "%f", &rinter);
832                                 break;
833 
834                         case DCM_IMGRESCALESLOPE:
835                                 fread(dimension, sizeof(char), dblock, fp_img);
836                                 dimension[dblock] = '\0';
837                                 sscanf(dimension, "%f", &rslope);
838                                 break;
839 
840                         case DCM_IMGROWS:
841                                 fread(&rows, sizeof(short), 1, fp_img);
842                                 short_swap((char *) &rows);
843                                 break;
844 
845                         case DCM_IMGCOLUMNS:
846                                 fread(&cols, sizeof(short), 1, fp_img);
847                                 short_swap((char *) &cols);
848                                 break;
849 
850                         case DCM_IMGPIXELREPRESENTATION:
851                                 fread(&sign, sizeof(short), 1, fp_img);
852                                 short_swap((char *) &sign);
853                                 break;
854 
855                         case DCM_IMGBITSSTORED:
856                                 fread(&sbits, sizeof(short), 1, fp_img);
857                                 short_swap((char *) &sbits);
858                                 break;
859 
860                         case DCM_IMGBITSALLOCATED:
861                                 fread(&abits, sizeof(short), 1, fp_img);
862                                 short_swap((char *) &abits);
863                                 break;
864 
865                         case DCM_DLMITEM:
866                         case DCM_DLMITEMDELIMITATIONITEM:
867                         case DCM_DLMSEQUENCEDELIMITATIONITEM:
868                                 break;
869 
870                         default:
871                                 for (i = 0; i < dblock; i++)
872                                         fread(&junk, sizeof(char), 1, fp_img);
873                                 break;
874                 }
875                 fread(&group, sizeof(short), 1, fp_img);
876                 short_swap((char *) &group);
877                 fread(&element, sizeof(short), 1, fp_img);
878                 short_swap((char *) &element);
879         }
880 
881         if (fp_img)
882         {
883                 /* 
884                  * move to start of image data depending on VR mode
885                  */
886                 fread(&dblock, sizeof(int), 1, fp_img);
887                 word_swap((char *) &dblock);
888                 dblock_vr_conv(&dblock, fp_img);
889 
890                 imregion.lx = 0;
891                 imregion.ux = cols;
892                 imregion.ly = 0;
893                 imregion.uy = rows;
894 
895                 if (sign == 0)
896                 {
897                         if (abits == 16 || abits == 12)
898                                 new_vtype = ushort_v;
899                         else
900                                 new_vtype = uchar_v;
901                 } else
902                 {
903                         if (abits == 16 || abits == 12)
904                                 new_vtype = short_v;
905                         else
906                                 new_vtype = char_v;
907                 }
908                 if (abits == 12)
909                         imregion.ux = 3.0 * imregion.ux / 4;
910                 cols = imregion.ux;
911 
912                 imrect = im_alloc(rows, cols, &imregion, (Vartype) new_vtype);
913 
914                 if (!fread_imrect_data(imrect, fp_img, pathname))
915                 {
916                         im_free(imrect);
917                         imrect = NULL;
918                         return (NULL);
919                 }
920                 if (abits == 12)
921                 {
922                         imrect2 = im_dicom_conv(imrect);
923                         im_free(imrect);
924                         imrect = imrect2;
925                 }
926                 if (imrect != NULL && get_swapping_ts())
927                 {
928                         im_endian_inplace(imrect);
929                 }
930 
931                 range = pow(2.0, sbits) - 1;
932 
933 printf("rinter: %4.2f  rslope: %4.2f\n", rinter, rslope);
934 
935                 imrect2 = im_cast(imrect, float_v);
936                 im_free(imrect);
937                 imrect = imrect2;
938                 if (rslope != 0.0 || rinter != 0.0)
939                 {
940                         for (j = imrect->region->ly; j < imrect->region->uy; j++)
941                                 for (k = imrect->region->lx; k < imrect->region->ux; k++)
942                                 {
943                                         pix = im_get_pixf(imrect, j, k);
944                                         pix = pix * rslope + rinter;
945                                         im_put_pixf(pix, imrect, j, k);
946                                 }
947                 }
948 
949                 iscale = vec3_alloc();
950                 iscale->el[0] = xsize;
951                 iscale->el[1] = ysize;
952                 iscale->el[2] = zsize;
953                 imrect->props = proplist_addifnp(imrect->props, iscale,
954                                                                                   VOXELS, vec3_free, false);
955         }
956 
957         return imrect;
958 }
959 
960 
961 
962 Imrect *dicom_read_image(const char *pathname)
963 {
964         FILE *fp_img;
965         Imrect *im;
966 
967         if ((fp_img = fopen_2(pathname, "rb")) == NULL)
968                 return (NULL);
969 
970         switch (dicom_preread_tests(fp_img))
971         {
972                 case PART10:
973                         im = dicom_read_multiformat_image(pathname, fp_img);
974                         fclose_2(fp_img, pathname);
975                         break;
976                 case NONPART10:
977                         im = dicom_read_multiformat_image(pathname, fp_img);
978                         fclose_2(fp_img, pathname);
979                         break;
980                 default:
981                         format("error: file failed pre checks\n");
982                         fclose_2(fp_img, pathname);
983                         im = NULL;
984                         break;
985         }
986 
987         return (im);
988 }
989 
990 
991 /*
992   recover information from dicom headers ref'd in seq using hdrextract_func
993   which should store info in image propslist.  hdrextract_func should be
994   passed fileptr and imrect and should return Bool success. 
995 */
996 Bool dicom_hdrinfo_extract(Sequence * seq, void (*hdrextract_func))
997 {
998         List *store = (List *) get_start_el();
999         Bool status = false, success = true;
1000         FILE *fp = NULL;
1001         char *pfname[1];
1002         char filename[512];
1003         char temp[512];
1004         int i;
1005 
1006         strcpy(filename, "");
1007 
1008         for (i = (seq->seq_start); i <= seq->seq_end; i++)
1009         {
1010                 parse_fname(seq->filename, temp, i);
1011                 sprintf(filename, "%s%s", temp, "");
1012                 *pfname = filename;
1013                 if (strchr(filename, '*') && !fname_resolve(NULL, filename, pfname))
1014                 {
1015                         error("dicom_hdrinfo_extract: no such file or filename", warning);
1016                         return (false);
1017                 }
1018 
1019                 if ((store == NULL) || (store->to == NULL))
1020                         continue;
1021                 if ((fp = fopen_2(*pfname, "rb")) == NULL)
1022                         continue;
1023                 status = ((Bool(*)())hdrextract_func) (fp, (Imrect *) (store->to));
1024                 fclose_2(fp, filename);
1025                 if (!status)
1026                         success = status;
1027                 store = store->next;
1028         }
1029 
1030         return (success);
1031 }
1032 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.