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

Linux Cross Reference
Tina6/tina-libs/tina/file/fileDicom_io.c

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

  1 /**********
  2  *
  3  * Copyright (c) 2003, Division of Imaging Science and Biomedical Engineering,
  4  * University of Manchester, UK.  All rights reserved.
  5  * 
  6  * Redistribution and use in source and binary forms, with or without modification, 
  7  * are permitted provided that the following conditions are met:
  8  * 
  9  *   . Redistributions of source code must retain the above copyright notice, 
 10  *     this list of conditions and the following disclaimer.
 11  *
 12  *   . Redistributions in binary form must reproduce the above copyright notice,
 13  *     this list of conditions and the following disclaimer in the documentation 
 14  *     and/or other materials provided with the distribution.
 15  * 
 16  *   . Neither the name of the University of Manchester nor the names of its
 17  *     contributors may be used to endorse or promote products derived from this 
 18  *     software without specific prior written permission.
 19  * 
 20  * 
 21  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
 22  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
 23  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
 24  * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
 25  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
 26  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
 27  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 28  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
 29  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
 30  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
 31  * POSSIBILITY OF SUCH DAMAGE.
 32  *
 33  **********
 34  *
 35  * Program :    TINA
 36  * File    :  $Source: /home/tina/cvs/tina-libs/tina/file/fileDicom_io.c,v $
 37  * Date    :  $Date: 2007/02/15 01:52:29 $
 38  * Version :  $Revision: 1.8 $
 39  * CVS Id  :  $Id: fileDicom_io.c,v 1.8 2007/02/15 01:52:29 paul Exp $
 40  *
 41  * Author  :  Legacy TINA
 42  *
 43  * Notes   :
 44  *
 45  * Functions to read/write DICOM image files
 46  *       These functions are a modified version of the
 47  *       ACR-NEMA ones
 48  *
 49  *********
 50 */
 51 
 52 #include "fileDicom_io.h"
 53 
 54 #if HAVE_CONFIG_H
 55   #include <config.h>
 56 #endif
 57 
 58 #include <math.h>
 59 #include <stdio.h>
 60 #include <string.h>
 61 #include <tina/sys/sysPro.h>
 62 #include <tina/sys/sysDef.h>
 63 #include <tina/image/imgPro.h>
 64 #include <tina/image/imgDef.h>
 65 #include <tina/math/mathPro.h>
 66 #include <tina/math/mathDef.h>
 67 #include <tina/file/file_NemaPro.h>
 68 #include <tina/file/file_DicomDef.h>
 69 #include <tina/file/file_UtilPro.h>
 70 
 71 
 72 #define VOXELS    450
 73 #define DYNSTIME  451
 74 #define TE_DATA   453
 75 #define PAT_DATA  454
 76 #define TR_DATA   455
 77 #define FLIP_ANGLE_DATA 456
 78 
 79 #define HEADERMAXCOUNT 100
 80 #define PART10    1
 81 #define NONPART10 0
 82 #define RUBBISH   -1
 83 
 84 #define I_LITTLE_ENDIAN "1.2.840.10008.1.2"
 85 #define E_LITTLE_ENDIAN "1.2.840.10008.1.2.1"
 86 #define E_BIG_ENDIAN    "1.2.840.10008.1.2.2"
 87 
 88 static int dicom_format = RUBBISH;
 89 
 90 static void set_dicom_format(int format)
 91 {
 92         dicom_format = format;
 93 }
 94 
 95 
 96 static int dblock_vr_conv(unsigned int *dblock, FILE * fp)
 97 {
 98         size_t readerr;
 99         unsigned int   dtemp;
100         unsigned char *pvr;
101 
102         if (dicom_format == RUBBISH)
103                 return (-1);
104 
105         pvr = (unsigned char *) dblock;
106 
107         /* 
108          * VR is a nested sequence (SQ)
109          */
110 
111         if ((pvr[0] == 'S' && pvr[1] == 'Q') || (pvr[3] == 'S' && pvr[2] == 'Q'))
112         {
113                 readerr = fread(dblock, sizeof(int), 1, fp);
114                 word_swap((char *) dblock);
115                 *dblock = 0;
116                 return (1);
117         }
118 
119         if ((pvr[0] == 'O' && pvr[1] == 'B') || (pvr[0] == 'O' && pvr[1] == 'W'))
120         {
121                 readerr = fread(dblock, sizeof(int), 1, fp);
122                 word_swap((char *) dblock);
123                 return (1);
124         }
125 
126         if ((pvr[3] == 'O' && pvr[2] == 'B') || (pvr[3] == 'O' && pvr[2] == 'W'))
127         {
128                 readerr = fread(dblock, sizeof(int), 1, fp);
129                 word_swap((char *) dblock);
130                 return (1);
131         }
132 
133         if (*dblock == 0xffffffff)
134         {
135                 *dblock = 0;
136                 return (1);
137         }
138 
139         if ((pvr[0] == 'A' && pvr[1] == 'E') ||
140                         (pvr[0] == 'A' && pvr[1] == 'S') ||
141                         (pvr[0] == 'A' && pvr[1] == 'T') ||
142                         (pvr[0] == 'C' && pvr[1] == 'S') ||
143                         (pvr[0] == 'D' && pvr[1] == 'A') ||
144                         (pvr[0] == 'D' && pvr[1] == 'S') ||
145                         (pvr[0] == 'D' && pvr[1] == 'T') ||
146                         (pvr[0] == 'F' && pvr[1] == 'L') ||
147                         (pvr[0] == 'F' && pvr[1] == 'D') ||
148                         (pvr[0] == 'I' && pvr[1] == 'S') ||
149                         (pvr[0] == 'L' && pvr[1] == 'O') ||
150                         (pvr[0] == 'L' && pvr[1] == 'T') ||
151                         (pvr[0] == 'P' && pvr[1] == 'N') ||
152                         (pvr[0] == 'S' && pvr[1] == 'H') ||
153                         (pvr[0] == 'S' && pvr[1] == 'L') ||
154                         (pvr[0] == 'S' && pvr[1] == 'S') ||
155                         (pvr[0] == 'S' && pvr[1] == 'T') ||
156                         (pvr[0] == 'T' && pvr[1] == 'M') ||
157                         (pvr[0] == 'U' && pvr[1] == 'I') ||
158                         (pvr[0] == 'U' && pvr[1] == 'L') || (pvr[0] == 'U' && pvr[1] == 'S'))
159         {
160                 /*
161                 (*dblock) &= 0xffff0000;
162                 (*dblock) >>= 16;
163                 */
164 
165 #ifdef LITTLE_ENDIAN_ARCHITECTURE
166                 dtemp = (unsigned int)(256*pvr[3] + pvr[2]);
167 #else
168                 dtemp = (unsigned int)(256*pvr[2] + pvr[3]);
169 #endif
170 
171                 (*dblock) = dtemp;
172                 return (1);
173         }
174 
175         if ((pvr[3] == 'A' && pvr[2] == 'E') ||
176                         (pvr[3] == 'A' && pvr[2] == 'S') ||
177                         (pvr[3] == 'A' && pvr[2] == 'T') ||
178                         (pvr[3] == 'C' && pvr[2] == 'S') ||
179                         (pvr[3] == 'D' && pvr[2] == 'A') ||
180                         (pvr[3] == 'D' && pvr[2] == 'S') ||
181                         (pvr[3] == 'D' && pvr[2] == 'T') ||
182                         (pvr[3] == 'F' && pvr[2] == 'L') ||
183                         (pvr[3] == 'F' && pvr[2] == 'D') ||
184                         (pvr[3] == 'I' && pvr[2] == 'S') ||
185                         (pvr[3] == 'L' && pvr[2] == 'O') ||
186                         (pvr[3] == 'L' && pvr[2] == 'T') ||
187                         (pvr[3] == 'P' && pvr[2] == 'N') ||
188                         (pvr[3] == 'S' && pvr[2] == 'H') ||
189                         (pvr[3] == 'S' && pvr[2] == 'L') ||
190                         (pvr[3] == 'S' && pvr[2] == 'S') ||
191                         (pvr[3] == 'S' && pvr[2] == 'T') ||
192                         (pvr[3] == 'T' && pvr[2] == 'M') ||
193                         (pvr[3] == 'U' && pvr[2] == 'I') ||
194                         (pvr[3] == 'U' && pvr[2] == 'L') || (pvr[3] == 'U' && pvr[2] == 'S'))
195         {
196                 /*
197                 (*dblock) &= 0x0000ffff;
198                 */
199 
200 #ifdef LITTLE_ENDIAN_ARCHITECTURE
201                 dtemp = (unsigned int)(256*pvr[1] + pvr[0]);
202 #else
203                 dtemp = (unsigned int)(256*pvr[0] + pvr[1]);
204 #endif
205 
206                 (*dblock) =  dtemp;
207                 return (1);
208         }
209 
210         /*
211                  must be implicit VR
212   */
213         return (-1);
214 }
215 
216 /*
217   finds endian UID in part10
218   leaves file ptr after last DCM_GROUPFILEMETA tag
219 */
220 static Bool dicom_part10_endian_swap(FILE * fp)
221 {
222         Bool success = false;
223         char endian[32];
224         unsigned char junk;
225         unsigned short group, element;
226         unsigned int i, dblock;
227         long int where = ftell(fp);
228 
229 /* part10 meta dicom are implicit little endian by default */
230 
231 #ifdef LITTLE_ENDIAN_ARCHITECTURE
232         Bool little = false;
233         Bool big = true;
234 
235         set_swapping_ts(0);
236 #else
237         Bool little = true;
238         Bool big = false;
239 
240         set_swapping_ts(1);
241 #endif
242 
243         fread(&group, sizeof(short), 1, fp);
244         short_swap((char *) &group);
245 
246         while (fp && group <= DCM_GROUPFILEMETA)
247         {
248                 fread(&element, sizeof(short), 1, fp);
249                 short_swap((char *) &element);
250                 fread(&dblock, sizeof(int), 1, fp);
251                 word_swap((char *) &dblock);
252                 dblock_vr_conv(&dblock, fp);
253 
254                 switch (DCM_MAKETAG(group, element))
255                 {
256                         case DCM_METATRANSFERSYNTAX:
257                                 fread(endian, sizeof(char), dblock, fp);
258                                 endian[dblock] = '\0';
259                                 success = true;
260                                 break;
261 
262                         case DCM_DLMITEM:
263                         case DCM_DLMITEMDELIMITATIONITEM:
264                         case DCM_DLMSEQUENCEDELIMITATIONITEM:
265                                 break;
266 
267                         default:
268                                 for (i = 0; i < dblock; i++)
269                                         fread(&junk, sizeof(char), 1, fp);
270                                 break;
271                 }
272 
273                 where = ftell(fp);
274                 fread(&group, sizeof(short), 1, fp);
275                 short_swap((char *) &group);
276         }
277 
278 /* put fptr back before last group read */
279 
280         fseek(fp, where, SEEK_SET);
281 
282         if (success && !strcmp(endian, E_BIG_ENDIAN))
283                 return (big);
284 
285         return (little);
286 }
287 
288 
289 /*
290   part 10 dicom has a preamble of 128 bytes followed by
291   "DICM", although preamble may be missing. 
292   sets endian.
293 */
294 static int dicom_preread_tests(FILE * fp)
295 {
296     size_t readerr;
297     short group, element, temp;
298     char td, ti, tc, tm;
299     int i, dblock, repeat, count;
300 
301     if (fp == NULL) return (RUBBISH);
302 
303     /* part10 DICM header */
304     fseek(fp, 128, SEEK_SET);
305     readerr = fread(&td, sizeof(char), 1, fp);
306     readerr = fread(&ti, sizeof(char), 1, fp);
307     readerr = fread(&tc, sizeof(char), 1, fp);
308     readerr = fread(&tm, sizeof(char), 1, fp);
309     if (td == 'D' && ti == 'I' && tc == 'C' && tm == 'M')
310     {
311         set_dicom_format(PART10);
312 
313         if (dicom_part10_endian_swap(fp))
314             set_swapping_ts(1);
315         else
316             set_swapping_ts(0);
317 
318         return (PART10);
319     }
320 
321     /* no pre-amble part10 DICM header */
322     rewind(fp);
323     readerr = fread(&td, sizeof(char), 1, fp);
324     readerr = fread(&ti, sizeof(char), 1, fp);
325     readerr = fread(&tc, sizeof(char), 1, fp);
326     readerr = fread(&tm, sizeof(char), 1, fp);
327     if (td == 'D' && ti == 'I' && tc == 'C' && tm == 'M')
328     {
329         set_dicom_format(PART10);
330 
331         if (dicom_part10_endian_swap(fp))
332             set_swapping_ts(1);
333         else
334             set_swapping_ts(0);
335 
336         return (PART10);
337     }
338 
339     repeat = 0;
340     set_swapping_ts(0);
341 
342     while (repeat < 2)
343     {
344         rewind(fp);
345         readerr = fread(&group, sizeof(short), 1, fp);
346         short_swap((char *) &group);
347         readerr = fread(&element, sizeof(short), 1, fp);
348         short_swap((char *) &element);
349         readerr = fread(&dblock, sizeof(int), 1, fp);
350         word_swap((char *) &dblock);
351 
352         count = 0;
353         while (group < DCM_GROUPIDENTIFYING && count < HEADERMAXCOUNT)
354         {
355             for (i = 0; i < dblock; i++)
356                 readerr = fread(&temp, sizeof(char), 1, fp);
357 
358             readerr = fread(&group, sizeof(short), 1, fp);
359             short_swap((char *) &group);
360             readerr = fread(&element, sizeof(short), 1, fp);
361             short_swap((char *) &element);
362             readerr = fread(&dblock, sizeof(int), 1, fp);
363             word_swap((char *) &dblock);
364             count++;
365         }
366 
367         /* standard non-part10 header start */
368         if (DCM_MAKETAG(group, element) == DCM_IDGROUPLENGTH && dblock == 4)
369         {
370             rewind(fp);
371             set_dicom_format(NONPART10);
372             return (NONPART10);
373         }
374 
375         /* inferior although common non-part10 header starts */
376         if ((DCM_MAKETAG(group, element) == DCM_IDLENGTHTOEND) ||
377             (DCM_MAKETAG(group, element) == DCM_IDIMAGETYPE) ||
378             (DCM_MAKETAG(group, element) == DCM_IDSPECIFICCHARACTER))
379         {
380             rewind(fp);
381             set_dicom_format(NONPART10);
382             return (NONPART10);
383         }
384 
385         set_swapping_ts(1);
386         repeat++;
387     }
388 
389     return (RUBBISH); /* failed */
390 }
391 
392 
393 Imrect *im_dicom_conv(Imrect * im)
394 {
395     Imrect *im2;
396     Imregion roi;
397     unsigned char *row1;
398     unsigned char *row2;
399     int i, j, k;
400     int lx, ux, ly, uy;
401 
402     if (im == NULL)
403         return (NULL);
404 
405     roi = *(im->region);
406     roi.ux = roi.ux * 4.0 / 3;
407     im2 = im_alloc(im->height, roi.ux, &roi, short_v);
408     lx = roi.lx;
409     ux = roi.ux;
410     ly = roi.ly;
411     uy = roi.uy;
412 
413     for (i = ly; i < uy; ++i)
414     {
415         row1 = IM_ROW(im, i);
416         row2 = IM_ROW(im2, i);
417         for (k = 2 * lx, j = 2 * lx; k < 2 * ux; j += 3)
418         {
419             row2[k++] = row1[j];
420             row2[k++] = row1[j + 1] % 16;
421             row2[k++] = row1[j + 1] / 16 + (row1[j + 2] % 16) * 16;
422             row2[k++] = row1[j + 2] / 16;
423         }
424     }
425 
426     return (im2);
427 }
428 
429 
430 Bool dicom_hdr_TE_extract()
431 {
432     float *TE_arr, *TR_arr, *FA_arr;
433     float te=0.0, tr=0.0, fa=0.0;
434     int type;
435     char temp[64];
436     unsigned char junk;
437     unsigned short group, element;
438     unsigned int dblock;
439 
440     Bool success = false;
441     FILE *fp = NULL;
442     char *pfname[1];
443     char filename[512];
444     char filename_temp[512];
445     int i, k;
446     Sequence *seq = (Sequence *)seq_get_current();
447     int end = get_end_frame(seq);
448 
449     strcpy(filename, "");
450 
451     /* These three vectors were originally allocated using fvector_alloc(seq->offset, end+1).   */
452     /* However, this causes a problem when the sequence is loaded from a non-zero offset (i.e.  */
453     /* a non-zero starting frame), since the proplist freeing function does not receive the     */
454     /* value of the offset, and so cannot set it back to zero before passing it to rfree. In    */
455     /* order to prevent the resultant crash, I have modified the allocation to have a zero      */
456     /* offset.  This should not cause any problems, as any function reading the vector would    */
457     /* have to take the offset into account anyway in order to read the correct element, so     */
458     /* the only change is a bit of wasted memory.  PAB 7/12/2011.                               */
459 
460     TE_arr = fvector_alloc(0, end+1);
461     TR_arr = fvector_alloc(0, end+1);
462     FA_arr = fvector_alloc(0, end+1);
463 
464     for (i = seq->offset; i <= end; i = i+seq->stride)
465     {
466         success=false;
467         parse_fname(seq->filename, filename_temp, i);
468         sprintf(filename, "%s%s", filename_temp, "");
469         *pfname = filename;
470         if (strchr(filename, '*') && !fname_resolve(NULL, filename, pfname))
471         {
472             error("dicom_hdrinfo_extract: no such file or filename", warning);
473             return (false);
474         }
475 
476         /*if ((store == NULL) || (store->to == NULL)) continue;*/
477 
478         if ((fp = fopen_2(*pfname, "rb")) == NULL)
479         {
480             error("file pointer equals NULL", warning);
481             return(false);
482         }
483 
484         if ((type = dicom_preread_tests(fp)) == RUBBISH)
485         {
486             format("error: file failed pre checks\n");
487             return (false);
488         }
489 
490         fread(&group, sizeof(short), 1, fp);
491         short_swap((char *) &group);
492 
493         while (fp /*&& !success*/ && group != DCM_GROUPPIXEL)
494         {
495             fread(&element, sizeof(short), 1, fp);
496             short_swap((char *) &element);
497             fread(&dblock, sizeof(int), 1, fp);
498             word_swap((char *) &dblock);
499             dblock_vr_conv(&dblock, fp);
500 
501             switch (DCM_MAKETAG(group, element))
502             {
503                 case DCM_ACQREPETITIONTIME:
504                     fread(temp, sizeof(char), dblock, fp);
505                     temp[dblock] = '\0';
506                     sscanf(temp, "%f", &tr);
507                     TR_arr[i]=tr;
508                 break;
509 
510                 case DCM_ACQECHOTIME:
511                     fread(temp, sizeof(char), dblock, fp);
512                     temp[dblock] = '\0';
513                     sscanf(temp, "%f", &te);
514                     TE_arr[i]=te;
515                 break;
516 
517                 case DCM_ACQFLIPANGLE:
518                     fread(temp, sizeof(char), dblock, fp);
519                     temp[dblock] = '\0';
520                     sscanf(temp, "%f", &fa);
521                     FA_arr[i]=fa;
522                 break;
523 
524                 case DCM_DLMITEM:
525                 case DCM_DLMITEMDELIMITATIONITEM:
526                 case DCM_DLMSEQUENCEDELIMITATIONITEM:
527                 break;
528 
529                 default:
530                     for (k = 0; k < dblock; k++)
531                     fread(&junk, sizeof(char), 1, fp);
532                 break;
533             }
534             fread(&group, sizeof(short), 1, fp);
535             short_swap((char *) &group);
536         }
537         fclose_2(fp, filename);
538     } 
539     seq->props = proplist_rm(seq->props, TR_DATA);
540     seq->props = proplist_addifnp(seq->props, (void *) TR_arr, TR_DATA, rfree, false);
541     seq->props = proplist_rm(seq->props, TE_DATA);
542     seq->props = proplist_addifnp(seq->props, (void *) TE_arr, TE_DATA, rfree, false);
543     seq->props = proplist_rm(seq->props, FLIP_ANGLE_DATA);
544     seq->props = proplist_addifnp(seq->props, (void *) FA_arr, FLIP_ANGLE_DATA, rfree, false);
545     success=true;
546     return (success); /*mjs 7/11/05 bit pointless as this just returns the last value of success, not really a measure of how well it all went. */
547 }
548 
549 
550 static int int_lut(char character)
551 {
552     /*mjs I would use atoi() but I can't seem to get it working */
553     if (character == 48) return 0;
554     if (character == 49) return 1;
555     if (character == 50) return 2;
556     if (character == 51) return 3;
557     if (character == 52) return 4;
558     if (character == 53) return 5;
559     if (character == 54) return 6;
560     if (character == 55) return 7;
561     if (character == 56) return 8;
562     if (character == 57) return 9;
563 
564     /* if (character == 32) */ /* If you get here, you already know the answer. PAB 9/1/2005*/ 
565     return 0;
566 }
567 
568 
569 static double get_dicom_time(char * temp)
570 {
571     /* cos DICOM is SHITE ... and so am i :)*/
572     double hour_1 =0.0, hour_2=0.0, min_1=0.0, min_2=0.0, time = 0.000000;
573     double secs_1=0.0, secs_2=0.0;
574     char h1=temp[0], h2=temp[1], m1=temp[2], m2=temp[3], s1=temp[4];
575 
576     hour_1 = (double)(int_lut(h1))*36000.0;
577     hour_2 = (double)(int_lut(h2))*3600.0;
578     min_1  = (double)(int_lut(m1))*600.0;
579     min_2  = (double)(int_lut(m2))*60.0;
580     secs_1 = (double)(int_lut(s1))*10.0;
581 
582     sscanf(&temp[5], "%lf", &secs_2);
583     time = secs_2 + secs_1 + min_2 + min_1 + hour_2 + hour_1;
584 
585     /*printf("%f %f %f %f %f %f\n", hour_1, hour_2, min_1, min_2, secs_1,secs_2);
586     printf("get_dicom_time: %f\n", time);*/
587     return(time);
588 }
589 
590 
591 /* Ok this function is quite frankly shafted */
592 Bool dicom_hdr_dynstimes_extract(Sequence *seq, int end)
593 {
594         double ser_time=0.0, im_time=0.0, time;
595         int type;
596         float *times;
597         char temp[64];
598         unsigned char junk;
599         unsigned short group, element;
600         unsigned int dblock;
601         Bool status = false, success = false;
602         FILE *fp = NULL;
603         char *pfname[1];
604         char filename[512];
605         char filename_temp[512];
606         int i, tend, k;
607         int j=0;
608 
609         strcpy(filename, "");
610 
611         tend = (end - seq->offset)+1;
612         
613         times = fvector_alloc(0, tend);
614         for (i = (seq->offset); i <= end; i = i+seq->stride)
615           {
616             status = false;
617             parse_fname(seq->filename, filename_temp, i);
618             sprintf(filename, "%s%s", filename_temp, "");
619             *pfname = filename;
620             if (strchr(filename, '*') && !fname_resolve(NULL, filename, pfname))
621               {
622                 error("dicom_hdrinfo_extract: no such file or filename", warning);
623                 return (false);
624               }
625 
626             /*if ((store == NULL) || (store->to == NULL))
627               continue;*/
628             if ((fp = fopen_2(*pfname, "rb")) == NULL)
629               {
630                 error("file pointer equals NULL", warning);
631                 return(false);
632               }
633             if ((type = dicom_preread_tests(fp)) == RUBBISH)
634               {
635                 format("error: file failed pre checks\n");
636                 return (false);
637               }
638 
639             fread(&group, sizeof(short), 1, fp);
640             short_swap((char *) &group);
641 
642             while (fp && !success && group != DCM_GROUPPIXEL)
643               {
644                 fread(&element, sizeof(short), 1, fp);
645                 short_swap((char *) &element);
646                 fread(&dblock, sizeof(int), 1, fp);
647                 word_swap((char *) &dblock);
648                 dblock_vr_conv(&dblock, fp);
649                 
650                 switch (DCM_MAKETAG(group, element))
651                   {
652                   case DCM_IDSERIESTIME:
653                     fread(temp, sizeof(char), dblock, fp);
654                     temp[dblock] = '\0';
655                     ser_time = get_dicom_time(temp);
656                     break;
657                     
658                   case DCM_IDIMAGETIME:
659                     fread(temp, sizeof(char), dblock, fp);
660                     temp[dblock] = '\0';
661                     im_time = get_dicom_time(temp);
662                     /*printf("image time    :%lf\n", im_time);*/
663                     success = true;
664                     break;
665                     
666                   case DCM_DLMITEM:
667                   case DCM_DLMITEMDELIMITATIONITEM:
668                   case DCM_DLMSEQUENCEDELIMITATIONITEM:
669                     break;
670                     
671                     
672                   default:
673                     for (k = 0; k < dblock; k++)
674                       fread(&junk, sizeof(char), 1, fp);
675                     break;
676                   }
677                 fread(&group, sizeof(short), 1, fp);
678                 short_swap((char *) &group);
679               }
680             
681             if (success)
682               {
683                 time = im_time - ser_time;
684                 times[j] =time;
685                 status = true;
686                 success = false;
687               }
688 
689             fclose_2(fp, filename);
690             j++;
691           }
692         if (status)
693           {
694             /* seq->props = proplist_rm(seq->props, DYNSTIME);
695                seq->props = proplist_addifnp(seq->props, (void *) times, DYNSTIME, dynstimes_free, false);*/
696 
697             seq->dim[3] = times[1];
698           }
699           
700         /* PAB 15/3/2011: Memory leak fixed: if we are not going to use the times vector, it should be freed */
701         fvector_free(times, 0);
702         
703         return (status);
704 }
705 
706 
707 Bool dicom_hdr_heartrate_extract(Sequence *seq, int end)
708 {
709   float *times=NULL, btime, phases, heartrate;
710   float beattime;
711   int type;
712   int psuccess = 0;
713   char temp[64];
714   unsigned char junk;
715   unsigned short group, element;
716   unsigned int dblock, tag1, tag2;
717   
718   Bool status = false, success = false;
719   FILE *fp = NULL;
720   char *pfname[1];
721   char filename[512];
722   char filename_temp[512];
723   int i, tend, k;
724   int j=0;
725 
726   strcpy(filename, "");
727 
728   /* mjs just temporarily, in futture extend with proper checks and stuff for 
729      all parameters */
730   tend = (end - seq->offset)+1;
731  
732   times = fvector_alloc(0, tend);
733  
734   for (i = (seq->offset); i <= end; i = i+seq->stride)
735     {
736       status = false;
737       parse_fname(seq->filename, filename_temp, i);
738       sprintf(filename, "%s%s", filename_temp, "");
739       *pfname = filename;
740       if (strchr(filename, '*') && !fname_resolve(NULL, filename, pfname))
741         {
742           error("dicom_hdrinfo_extract: no such file or filename", warning);
743           return (false);
744         }
745       
746       /*if ((store == NULL) || (store->to == NULL))
747         continue;*/
748       if ((fp = fopen_2(*pfname, "rb")) == NULL)
749         {       
750           error("file pointer equals NULL", warning);
751           return(false);
752         }
753       if ((type = dicom_preread_tests(fp)) == RUBBISH)
754         {
755           format("error: file failed pre checks\n");
756           return (false);
757         }
758       
759       fread(&group, sizeof(short), 1, fp);
760       short_swap((char *) &group);
761       
762       tag1 = 0x200117; /* ie, private group 2001,xx17 */
763       
764       while (fp && !success && group != DCM_GROUPPIXEL)
765         {
766           fread(&element, sizeof(short), 1, fp);
767           short_swap((char *) &element);
768           fread(&dblock, sizeof(int), 1, fp);
769           word_swap((char *) &dblock);
770           dblock_vr_conv(&dblock, fp);
771           
772           if (group == 0x2001)
773             {
774               format("0x2001 found");
775             }
776           tag2 = (group << 8) | (element & 0xFF); 
777           
778           if (tag2 == tag1) /* looks like it's 1017 anyway*/
779             {
780               fread(temp, sizeof(char), dblock, fp);
781               temp[dblock] = '\0';
782               sscanf(temp, "%f", &phases);
783               if (psuccess == 2)
784                 success = true;
785               else
786                 psuccess++;
787             }
788           else
789             {
790               
791               switch (DCM_MAKETAG(group, element))
792                 {/*
793                   case DCM_MAKETAG(0x0019, 0x1022):
794                   fread(temp, sizeof(char), dblock, fp);
795                   temp[dblock] = '\0';
796                   sscanf(temp, "%f", &btime);
797                   if (psuccess == 2)
798                     success = true;
799                   else
800                     psuccess++;
801                     break;*/
802                   /* mjs change no of phases tag - need to check*/
803                                 /* if 0x10 is right */
804                   
805                                 /* mjs change heart rate tag*/
806                 case DCM_MAKETAG(0x0018, 0x1088):
807                   fread(temp, sizeof(char), dblock, fp);
808                   temp[dblock] = '\0';
809                   sscanf(temp, "%f", &heartrate);
810                   /* if (psuccess == 1)*/
811                   success = true;
812                   /*else
813                     psuccess++;*/
814                   break;
815                   
816                 case DCM_DLMITEM:
817                 case DCM_DLMITEMDELIMITATIONITEM:
818                 case DCM_DLMSEQUENCEDELIMITATIONITEM:
819                   break;
820                   
821                   
822                 default:
823                   for (k = 0; k < dblock; k++)
824                     fread(&junk, sizeof(char), 1, fp);
825                   break;
826                 }
827             }
828           fread(&group, sizeof(short), 1, fp);
829           short_swap((char *) &group);
830         }
831       phases = 15;
832       if (phases == 1)
833         {
834           success = false;
835         }
836       if (success)
837         {
838           btime = 0.0;
839           beattime = 1.0 / (heartrate / 60);
840           times[j] = btime + (float) j *(beattime / (phases - 1));
841           success = false;
842           status = true;
843         }
844          
845       fclose_2(fp, filename);
846       j++;
847     }
848   if (status)
849     {
850       /*    seq->props = proplist_rm(seq->props, DYNSTIME);
851             seq->props = proplist_addifnp(seq->props, (void *) times, DYNSTIME, dynstimes_free, false);  */
852       seq->dim[3] = times[1];
853     }
854   
855   return (status);
856   
857 }
858 
859 
860 Bool dicom_hdr_patientdetails_extract()
861 {/* again, can assume that the details are the same for all images */
862   Bool success = false;
863   int type;
864   char *details = NULL;
865   unsigned char junk;
866   unsigned short group, element;
867   unsigned int dblock;
868   FILE *fp = NULL;
869   char *pfname[1];
870   char filename[512];
871   char filename_temp[512];
872   int  k;
873   Sequence *seq = (Sequence *)seq_get_current();
874 
875   strcpy(filename, "");
876    
877   parse_fname(seq->filename, filename_temp, seq->offset);
878   sprintf(filename, "%s%s", filename_temp, "");
879   *pfname = filename;
880   if (strchr(filename, '*') && !fname_resolve(NULL, filename, pfname))
881     {
882       error("dicom_hdrinfo_extract: no such file or filename", warning);
883       return (false);
884     }
885   
886   if ((fp = fopen_2(*pfname, "rb")) == NULL)
887     {
888       error("file pointer equals NULL", warning);
889       return(false);
890     }
891   if ((type = dicom_preread_tests(fp)) == RUBBISH)
892     {
893       format("error: file failed pre checks\n");
894       return (false);
895     }
896   
897   fread(&group, sizeof(short), 1, fp);
898   short_swap((char *) &group);
899   
900   details = (char *) ralloc(256*sizeof(char));
901   while (fp && !success && group != DCM_GROUPPIXEL)
902     {
903       fread(&element, sizeof(short), 1, fp);
904       short_swap((char *) &element);
905       fread(&dblock, sizeof(int), 1, fp);
906       word_swap((char *) &dblock);
907       dblock_vr_conv(&dblock, fp);
908       
909       switch (DCM_MAKETAG(group, element))
910         {
911         case DCM_MAKETAG(0x0008, 0x0020):
912           fread(details, sizeof(char), dblock, fp);
913           details[dblock] = '\0';
914           success = true;
915           break;
916           
917         case DCM_DLMITEM:
918         case DCM_DLMITEMDELIMITATIONITEM:
919         case DCM_DLMSEQUENCEDELIMITATIONITEM:
920           break;
921           
922           
923         default:
924           for (k = 0; k < dblock; k++)
925             fread(&junk, sizeof(char), 1, fp);
926           break;
927         }
928       fread(&group, sizeof(short), 1, fp);
929       short_swap((char *) &group);
930     
931     }
932   if (success != true)
933     cvector_free(details, 0);
934   else
935     {
936       seq->props = proplist_rm(seq->props, PAT_DATA);
937       seq->props = proplist_addifnp(seq->props, (void *) details, PAT_DATA, rfree, false);
938     }
939   
940   fclose_2(fp, filename);
941   return (success);
942   
943 }
944 
945 
946 Bool dicom_hdr_voxelscale_extract()
947 {/*logically only need to look at the first image. */
948     float xsize, ysize, zsize, zzsize;
949     int type;
950     int xysuccess = 0, zsuccess = 0, zzsuccess = 0;
951     char dimension[64];
952     unsigned char junk;
953     unsigned short group, element;
954     unsigned int dblock;
955     Bool success = false;
956     FILE *fp = NULL;
957     char *pfname[1];
958     char filename[512];
959     char filename_temp[512];
960     int k;
961     Sequence * seq = (Sequence *)seq_get_current();
962 
963     strcpy(filename, "");
964 
965     parse_fname(seq->filename, filename_temp, seq->offset);
966     sprintf(filename, "%s%s", filename_temp, "");
967     *pfname = filename;
968 
969     if (strchr(filename, '*') && !fname_resolve(NULL, filename, pfname))
970     {
971         error("dicom_hdrinfo_extract: no such file or filename", warning);
972         return (false);
973     }
974 
975     if ((fp = fopen_2(*pfname, "rb")) == NULL)
976     {
977         error("file pointer equals NULL", warning);
978         return(false);
979     }
980     xsize = ysize = zsize = 1.0;
981 
982     if ((type = dicom_preread_tests(fp)) == RUBBISH)
983     {
984         format("error: file failed pre checks\n");
985         return (false);
986     }
987 
988     fread(&group, sizeof(short), 1, fp);
989     short_swap((char *) &group);
990 
991     while (fp && group != DCM_GROUPPIXEL)
992     {
993         fread(&element, sizeof(short), 1, fp);
994         short_swap((char *) &element);
995         fread(&dblock, sizeof(int), 1, fp);
996         word_swap((char *) &dblock);
997         dblock_vr_conv(&dblock, fp);
998 
999         switch (DCM_MAKETAG(group, element))
1000         {
1001             case DCM_ACQSLICESPACING:
1002                 fread(dimension, sizeof(char), dblock, fp);
1003                 dimension[dblock] = '\0';
1004                 sscanf(dimension, "%f", &zsize);
1005                 zsuccess++;
1006             break;
1007 
1008             case DCM_ACQSLICETHICKNESS:
1009                 fread(dimension, sizeof(char), dblock, fp);
1010                 dimension[dblock] = '\0';
1011                 sscanf(dimension, "%f", &zzsize);
1012                 zzsuccess++;
1013             break;
1014 
1015             case DCM_IMGPIXELSPACING:
1016                 fread(dimension, sizeof(char), dblock, fp);
1017                 dimension[dblock] = '\0';
1018                 sscanf(dimension, "%f\\%f", &xsize, &ysize);
1019                 xysuccess++;
1020             break;
1021 
1022             case DCM_DLMITEM:
1023             case DCM_DLMITEMDELIMITATIONITEM:
1024             case DCM_DLMSEQUENCEDELIMITATIONITEM:
1025             break;
1026 
1027             default:
1028                 for (k = 0; k < dblock; k++)
1029                 fread(&junk, sizeof(char), 1, fp);
1030             break;
1031         }
1032         fread(&group, sizeof(short), 1, fp);
1033         short_swap((char *) &group);
1034     }
1035 
1036     if (xysuccess==1)
1037     {
1038         seq->dim[0] = xsize;
1039         seq->dim[1] = ysize;
1040 
1041         if(zsuccess==1)
1042         {
1043             seq->dim[2] = zsize;
1044             success = true;
1045         }
1046         else if(zzsuccess==1)
1047         {
1048             seq->dim[2] = zzsize;
1049             success = true;
1050         }
1051     }
1052 
1053     fclose_2(fp, filename);
1054     return success;
1055 }
1056 
1057 
1058 static Imrect *dicom_read_multiformat_image(const char *pathname, FILE * fp_img, int scale_factor_flag)
1059 {
1060         int get_swapping_ts();
1061         Imrect *imrect = NULL, *imrect2 = NULL;
1062         Imregion imregion;
1063         Vartype new_vtype;
1064         float range;
1065         float rinter, rslope;
1066         float sinter, sslope; /*mjs added 25/11/05 to give functionality for loading scale slope etc */
1067         float pix;
1068         int i, j, k;
1069         unsigned char junk;
1070         char dimension[64];
1071         unsigned short abits;
1072         unsigned short sbits;
1073         unsigned short sign;
1074         unsigned short rows;
1075         unsigned short cols;
1076         unsigned short group;
1077         unsigned short element;
1078         unsigned int dblock;
1079         float xsize, ysize, zsize;
1080         Sequence *seq = (Sequence *) seq_get_current();
1081 
1082         rinter = 0.0;
1083         rslope = 1.0;
1084         xsize = ysize = zsize = 1.0;
1085 
1086         fread(&group, sizeof(short), 1, fp_img);
1087         short_swap((char *) &group);
1088         fread(&element, sizeof(short), 1, fp_img);
1089         short_swap((char *) &element);
1090 
1091         while (fp_img && DCM_MAKETAG(group, element) != DCM_PXLPIXELDATA)
1092         {
1093                 fread(&dblock, sizeof(int), 1, fp_img);
1094                 word_swap((char *) &dblock);
1095                 dblock_vr_conv(&dblock, fp_img);
1096 
1097                 if (DCM_MAKETAG(group, element) == DCM_PXLPIXELDATA)
1098                         break;
1099 
1100                 switch (DCM_MAKETAG(group, element))
1101                 {
1102                         case DCM_ACQSLICESPACING:
1103                                 fread(dimension, sizeof(char), dblock, fp_img);
1104                                 dimension[dblock] = '\0';
1105                                 sscanf(dimension, "%f", &zsize);
1106                                 break;
1107 
1108                         case DCM_IMGPIXELSPACING:
1109                                 fread(dimension, sizeof(char), dblock, fp_img);
1110                                 dimension[dblock] = '\0';
1111                                 sscanf(dimension, "%f\\%f", &xsize, &ysize);
1112                                 break;
1113 
1114                         case DCM_IMGRESCALEINTERCEPT:
1115                                 fread(dimension, sizeof(char), dblock, fp_img);
1116                                 dimension[dblock] = '\0';
1117                                 sscanf(dimension, "%f", &rinter);
1118                                 break;
1119 
1120                                 /*case DCM_IMGRESCALESLOPE:*/
1121                 case DCM_MAKETAG(0x0028,0x1053):        
1122                                 fread(dimension, sizeof(char), dblock, fp_img);
1123                                 dimension[dblock] = '\0';
1124                                 sscanf(dimension, "%f", &rslope);
1125                                 break;
1126                                 /*mjs 25/11/05 added following two tags in order to load in scale slope and intercept. note entertainingly they are Float VR's and not digital strings.*/
1127                         case DCM_MAKETAG(0x2005,0x100d):
1128                           fread(&sinter, sizeof(float), 1, fp_img);
1129                           break;
1130                         case DCM_MAKETAG(0x2005,0x100e):
1131                     
1132                           fread(&sslope, sizeof(float), 1, fp_img);
1133                           break;
1134                                 
1135                         case DCM_IMGROWS:
1136                                 fread(&rows, sizeof(short), 1, fp_img);
1137                                 short_swap((char *) &rows);
1138                                 break;
1139 
1140                         case DCM_IMGCOLUMNS:
1141                                 fread(&cols, sizeof(short), 1, fp_img);
1142                                 short_swap((char *) &cols);
1143                                 break;
1144 
1145                         case DCM_IMGPIXELREPRESENTATION:
1146                                 fread(&sign, sizeof(short), 1, fp_img);
1147                                 short_swap((char *) &sign);
1148                                 break;
1149 
1150                         case DCM_IMGBITSSTORED:
1151                                 fread(&sbits, sizeof(short), 1, fp_img);
1152                                 short_swap((char *) &sbits);
1153                                 break;
1154 
1155                         case DCM_IMGBITSALLOCATED:
1156                                 fread(&abits, sizeof(short), 1, fp_img);
1157                                 short_swap((char *) &abits);
1158                                 break;
1159 
1160                         case DCM_DLMITEM:
1161                         case DCM_DLMITEMDELIMITATIONITEM:
1162                         case DCM_DLMSEQUENCEDELIMITATIONITEM:
1163                                 break;
1164 
1165                         default:
1166                                 for (i = 0; i < dblock; i++)
1167                                         fread(&junk, sizeof(char), 1, fp_img);
1168                                 break;
1169                 }
1170                 fread(&group, sizeof(short), 1, fp_img);
1171                 short_swap((char *) &group);
1172                 fread(&element, sizeof(short), 1, fp_img);
1173                 short_swap((char *) &element);
1174         }
1175 
1176         if (fp_img)
1177         {
1178                 /* 
1179                  * move to start of image data depending on VR mode
1180                  */
1181                 fread(&dblock, sizeof(int), 1, fp_img);
1182                 word_swap((char *) &dblock);
1183                 dblock_vr_conv(&dblock, fp_img);
1184 
1185                 imregion.lx = 0;
1186                 imregion.ux = cols;
1187                 imregion.ly = 0;
1188                 imregion.uy = rows;
1189 
1190                 if (sign == 0)
1191                 {
1192                         if (abits == 16 || abits == 12)
1193                                 new_vtype = ushort_v;
1194                         else
1195                                 new_vtype = uchar_v;
1196                 } else
1197                 {
1198                         if (abits == 16 || abits == 12)
1199                                 new_vtype = short_v;
1200                         else
1201                                 new_vtype = char_v;
1202                 }
1203                 if (abits == 12)
1204                         imregion.ux = 3.0 * imregion.ux / 4;
1205                 cols = imregion.ux;
1206 
1207                 imrect = im_alloc(rows, cols, &imregion, (Vartype) new_vtype);
1208 
1209                 if (!fread_imrect_data(imrect, fp_img, pathname))
1210                 {
1211                         im_free(imrect);
1212                         imrect = NULL;
1213                         return (NULL);
1214                 }
1215                 if (abits == 12)
1216                 {
1217                         imrect2 = im_dicom_conv(imrect);
1218                         im_free(imrect);
1219                         imrect = imrect2;
1220                 }
1221                 if (imrect != NULL && get_swapping_ts())
1222                 {
1223                         im_endian_inplace(imrect);
1224                 }
1225 
1226                 range = pow(2.0, sbits) - 1;
1227 
1228 
1229                 imrect2 = im_cast(imrect, float_v);
1230                 im_free(imrect);
1231                 imrect = imrect2;
1232                 if (rslope != 0.0 || rinter != 0.0)
1233                 {
1234                         for (j = imrect->region->ly; j < imrect->region->uy; j++)
1235                                 for (k = imrect->region->lx; k < imrect->region->ux; k++)
1236                                 {
1237                                         pix = im_get_pixf(imrect, j, k);
1238                                         if (scale_factor_flag == 0)
1239                                           {
1240                                             pix = pix * rslope + rinter;
1241                                           }
1242                                         if (scale_factor_flag == 1)
1243                                           {
1244                                             pix = (pix-sinter)/sslope;
1245                                           }
1246                                         im_put_pixf(pix, imrect, j, k);
1247                                 }
1248                 }
1249         
1250                   if (seq != NULL)
1251                   {
1252                           seq->dim[0] = xsize;
1253                           seq->dim[1] = ysize;
1254                           seq->dim[2] = zsize;
1255                   }
1256                 
1257         }
1258 
1259         return imrect;
1260 }
1261 
1262 
1263 
1264 Imrect *dicom_read_image(char *pathname, int scale_factor_flag)
1265 {
1266         FILE *fp_img;
1267         Imrect *im;
1268 
1269         if ((fp_img = fopen_2(pathname, "rb")) == NULL)
1270                 return (NULL);
1271 
1272         switch (dicom_preread_tests(fp_img))
1273         {
1274                 case PART10:
1275                         im = dicom_read_multiformat_image(pathname, fp_img, scale_factor_flag);
1276                         fclose_2(fp_img, pathname);
1277                         break;
1278                 case NONPART10:
1279                         im = dicom_read_multiformat_image(pathname, fp_img, scale_factor_flag);
1280                         fclose_2(fp_img, pathname);
1281                         break;
1282                 default:
1283                         format("error: file failed pre checks\n");
1284                         fclose_2(fp_img, pathname);
1285                         im = NULL;
1286                         break;
1287         }
1288 
1289         return (im);
1290 }
1291 
1292 
1293 /*
1294   recover information from dicom headers ref'd in seq using hdrextract_func
1295   which should store info in image propslist.  hdrextract_func should be
1296   passed fileptr and imrect and should return Bool success. 
1297 */
1298 /*
1299 Bool dicom_hdrinfo_extract(Sequence * seq, void (*hdrextract_func))
1300 {
1301         List *store = (List *) get_seq_start_el(seq);
1302         Bool status = false, success = true;
1303         FILE *fp = NULL;
1304         char *pfname[1];
1305         char filename[512];
1306         char temp[512];
1307         int i;
1308         int end = get_end_frame(seq);
1309 
1310         strcpy(filename, "");
1311 
1312         for (i = (seq->offset); i <= end; i = i+seq->stride)
1313         {
1314                 parse_fname(seq->filename, temp, i);
1315                 sprintf(filename, "%s%s", temp, "");
1316                 *pfname = filename;
1317                 if (strchr(filename, '*') && !fname_resolve(NULL, filename, pfname))
1318                 {
1319                         error("dicom_hdrinfo_extract: no such file or filename", warning);
1320                         return (false);
1321                 }
1322 
1323                 if ((store == NULL) || (store->to == NULL))
1324                         continue;
1325                 if ((fp = fopen_2(*pfname, "rb")) == NULL)
1326                         continue;
1327                 status = ((Bool(*)())hdrextract_func) (fp, (Imrect *) (store->to));
1328                 fclose_2(fp, filename);
1329                 if (!status)
1330                         success = status;
1331                 store = store->next;
1332         }
1333 
1334         return (success);
1335 }
1336 */
1337 

~ [ 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.