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

Linux Cross Reference
Tina5/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)
302                 return (RUBBISH);
303 
304 /* part10 DICM header */
305         fseek(fp, 128, SEEK_SET);
306         readerr = fread(&td, sizeof(char), 1, fp);
307         readerr = fread(&ti, sizeof(char), 1, fp);
308         readerr = fread(&tc, sizeof(char), 1, fp);
309         readerr = fread(&tm, sizeof(char), 1, fp);
310         if (td == 'D' && ti == 'I' && tc == 'C' && tm == 'M')
311         {
312                 set_dicom_format(PART10);
313 
314                 if (dicom_part10_endian_swap(fp))
315                         set_swapping_ts(1);
316                 else
317                         set_swapping_ts(0);
318 
319                 return (PART10);
320         }
321 
322 /* no pre-amble part10 DICM header */
323         rewind(fp);
324         readerr = fread(&td, sizeof(char), 1, fp);
325         readerr = fread(&ti, sizeof(char), 1, fp);
326         readerr = fread(&tc, sizeof(char), 1, fp);
327         readerr = fread(&tm, sizeof(char), 1, fp);
328         if (td == 'D' && ti == 'I' && tc == 'C' && tm == 'M')
329         {
330                 set_dicom_format(PART10);
331 
332                 if (dicom_part10_endian_swap(fp))
333                         set_swapping_ts(1);
334                 else
335                         set_swapping_ts(0);
336 
337                 return (PART10);
338         }
339 
340         repeat = 0;
341         set_swapping_ts(0);
342 
343         while (repeat < 2)
344         {
345                 rewind(fp);
346                 readerr = fread(&group, sizeof(short), 1, fp);
347                 short_swap((char *) &group);
348                 readerr = fread(&element, sizeof(short), 1, fp);
349                 short_swap((char *) &element);
350                 readerr = fread(&dblock, sizeof(int), 1, fp);
351                 word_swap((char *) &dblock);
352 
353                 count = 0;
354                 while (group < DCM_GROUPIDENTIFYING && count < HEADERMAXCOUNT)
355                 {
356                         for (i = 0; i < dblock; i++)
357                                 readerr = fread(&temp, sizeof(char), 1, fp);
358 
359                         readerr = fread(&group, sizeof(short), 1, fp);
360                         short_swap((char *) &group);
361                         readerr = fread(&element, sizeof(short), 1, fp);
362                         short_swap((char *) &element);
363                         readerr = fread(&dblock, sizeof(int), 1, fp);
364                         word_swap((char *) &dblock);
365                         count++;
366                 }
367 
368 /* standard non-part10 header start */
369                 if (DCM_MAKETAG(group, element) == DCM_IDGROUPLENGTH && dblock == 4)
370                 {
371                         rewind(fp);
372                         set_dicom_format(NONPART10);
373                         return (NONPART10);
374                 }
375 
376 /* inferior although common non-part10 header starts */
377          if ((DCM_MAKETAG(group, element) == DCM_IDLENGTHTOEND) ||
378              (DCM_MAKETAG(group, element) == DCM_IDIMAGETYPE) ||
379              (DCM_MAKETAG(group, element) == DCM_IDSPECIFICCHARACTER))
380                 {
381                         rewind(fp);
382                         set_dicom_format(NONPART10);
383                         return (NONPART10);
384                 }
385 
386                 set_swapping_ts(1);
387                 repeat++;
388         }
389 
390         return (RUBBISH);                                                       /* failed */
391 }
392 
393 
394 Imrect *im_dicom_conv(Imrect * im)
395 {
396         Imrect *im2;
397         Imregion roi;
398         unsigned char *row1;
399         unsigned char *row2;
400         int i, j, k;
401         int lx, ux, ly, uy;
402 
403         if (im == NULL)
404                 return (NULL);
405 
406         roi = *(im->region);
407         roi.ux = roi.ux * 4.0 / 3;
408         im2 = im_alloc(im->height, roi.ux, &roi, short_v);
409         lx = roi.lx;
410         ux = roi.ux;
411         ly = roi.ly;
412         uy = roi.uy;
413 
414         for (i = ly; i < uy; ++i)
415         {
416                 row1 = IM_ROW(im, i);
417                 row2 = IM_ROW(im2, i);
418                 for (k = 2 * lx, j = 2 * lx; k < 2 * ux; j += 3)
419                 {
420                         row2[k++] = row1[j];
421                         row2[k++] = row1[j + 1] % 16;
422                         row2[k++] = row1[j + 1] / 16 + (row1[j + 2] % 16) * 16;
423                         row2[k++] = row1[j + 2] / 16;
424                 }
425         }
426 
427         return (im2);
428 }
429 
430 
431 Bool dicom_hdr_TE_extract()
432 {
433         float *TE_arr, *TR_arr, *FA_arr;
434         float te=0.0, tr=0.0, fa=0.0;
435         int type;
436         char temp[64];
437         unsigned char junk;
438         unsigned short group, element;
439         unsigned int dblock;
440 
441         Bool success = false;
442         FILE *fp = NULL;
443         char *pfname[1];
444         char filename[512];
445         char filename_temp[512];
446         int i, k;
447         Sequence *seq = (Sequence *)seq_get_current();
448         int end = get_end_frame(seq);
449  
450         strcpy(filename, "");
451 
452         TE_arr = fvector_alloc(seq->offset, end+1);
453         TR_arr = fvector_alloc(seq->offset, end+1);
454         FA_arr = fvector_alloc(seq->offset, end+1);
455 
456         for (i = seq->offset; i <= end; i = i+seq->stride)
457         {
458         success=false;
459         parse_fname(seq->filename, filename_temp, i);
460         sprintf(filename, "%s%s", filename_temp, "");
461         *pfname = filename;
462         if (strchr(filename, '*') && !fname_resolve(NULL, filename, pfname))
463           {
464             error("dicom_hdrinfo_extract: no such file or filename", warning);
465             return (false);
466           }
467         
468         /*if ((store == NULL) || (store->to == NULL))
469           continue;*/
470         if ((fp = fopen_2(*pfname, "rb")) == NULL)
471           {
472             error("file pointer equals NULL", warning);
473             return(false);
474           }         
475         
476         if ((type = dicom_preread_tests(fp)) == RUBBISH)
477           {
478             format("error: file failed pre checks\n");
479             return (false);
480           }
481         
482         /*TE = (float *) ralloc(sizeof(float));
483         TR = (float *) ralloc(sizeof(float));
484         FA = (float *) ralloc(sizeof(float));*/
485 
486         fread(&group, sizeof(short), 1, fp);
487         short_swap((char *) &group);
488         
489         while (fp /*&& !success*/ && group != DCM_GROUPPIXEL)
490           {
491             fread(&element, sizeof(short), 1, fp);
492             short_swap((char *) &element);
493             fread(&dblock, sizeof(int), 1, fp);
494             word_swap((char *) &dblock);
495             dblock_vr_conv(&dblock, fp);
496             
497             switch (DCM_MAKETAG(group, element))
498               {         
499               case DCM_ACQREPETITIONTIME:
500                 fread(temp, sizeof(char), dblock, fp);
501                 temp[dblock] = '\0';
502                 sscanf(temp, "%f", &tr);
503                 TR_arr[i]=tr;
504                 break;
505 
506               case DCM_ACQECHOTIME:
507                 fread(temp, sizeof(char), dblock, fp);
508                 temp[dblock] = '\0';
509                 sscanf(temp, "%f", &te);
510                 TE_arr[i]=te;
511                 break;
512 
513               case DCM_ACQFLIPANGLE:
514                 fread(temp, sizeof(char), dblock, fp);
515                 temp[dblock] = '\0';
516                 sscanf(temp, "%f", &fa);
517                 FA_arr[i]=fa;
518                 break;
519 
520 
521               case DCM_DLMITEM:
522               case DCM_DLMITEMDELIMITATIONITEM:
523               case DCM_DLMSEQUENCEDELIMITATIONITEM:
524                 break;
525                 
526               default:
527                 for (k = 0; k < dblock; k++)
528                   fread(&junk, sizeof(char), 1, fp);
529                 break;
530               }
531             fread(&group, sizeof(short), 1, fp);
532             short_swap((char *) &group);
533           }           
534         fclose_2(fp, filename);
535 
536         } 
537         seq->props = proplist_rm(seq->props, TR_DATA);
538         seq->props = proplist_addifnp(seq->props, (void *) TR_arr, TR_DATA, rfree, false);
539         seq->props = proplist_rm(seq->props, TE_DATA);
540         seq->props = proplist_addifnp(seq->props, (void *) TE_arr, TE_DATA, rfree, false);
541         seq->props = proplist_rm(seq->props, FLIP_ANGLE_DATA);
542         seq->props = proplist_addifnp(seq->props, (void *) FA_arr, FLIP_ANGLE_DATA, rfree, false);
543         success=true;
544         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. */
545         
546 }
547 
548 
549 static int int_lut(char character)
550 {
551   /*mjs I would use atoi() but I can't seem to get it working */
552   if (character == 48)
553     return 0;
554   if (character == 49)
555     return 1;
556   if (character == 50)
557     return 2;
558   if (character == 51)
559     return 3;
560   if (character == 52)
561     return 4;
562   if (character == 53)
563     return 5;
564   if (character == 54)
565     return 6;
566   if (character == 55)
567     return 7;
568   if (character == 56)
569     return 8;
570   if (character == 57)
571     return 9;
572   /* if (character == 32) */ /* If you get here, you already know the answer. PAB 9/1/2005*/ 
573     return 0;
574 }
575 
576 
577 static double get_dicom_time(char * temp)
578 {
579   /* cos DICOM is SHITE ... and so am i :)*/
580   double hour_1 =0.0, hour_2=0.0, min_1=0.0, min_2=0.0, time = 0.000000;
581   double secs_1=0.0, secs_2=0.0;
582   char h1=temp[0], h2=temp[1], m1=temp[2], m2=temp[3], s1=temp[4];
583 
584   hour_1 = (double)(int_lut(h1))*36000.0;
585   hour_2 = (double)(int_lut(h2))*3600.0;
586   min_1  = (double)(int_lut(m1))*600.0;
587   min_2  = (double)(int_lut(m2))*60.0;
588   secs_1 = (double)(int_lut(s1))*10.0;
589   
590   sscanf(&temp[5], "%lf", &secs_2);
591   time = secs_2 + secs_1 + min_2 + min_1 + hour_2 + hour_1;
592   
593   /*printf("%f %f %f %f %f %f\n", hour_1, hour_2, min_1, min_2, secs_1,secs_2);
594     printf("get_dicom_time: %f\n", time);*/
595   return(time);
596 }
597 
598 /* Ok this function is quite frankly shafted */
599 Bool dicom_hdr_dynstimes_extract(Sequence *seq, int end)
600 {
601         double ser_time=0.0, im_time=0.0, time;
602         int type;
603         float *times;
604         char temp[64];
605         unsigned char junk;
606         unsigned short group, element;
607         unsigned int dblock;
608         Bool status = false, success = false;
609         FILE *fp = NULL;
610         char *pfname[1];
611         char filename[512];
612         char filename_temp[512];
613         int i, tend, k;
614         int j=0;
615 
616         strcpy(filename, "");
617 
618         tend = (end - seq->offset)+1;
619         
620         times = fvector_alloc(0, tend);
621         for (i = (seq->offset); i <= end; i = i+seq->stride)
622           {
623             status = false;
624             parse_fname(seq->filename, filename_temp, i);
625             sprintf(filename, "%s%s", filename_temp, "");
626             *pfname = filename;
627             if (strchr(filename, '*') && !fname_resolve(NULL, filename, pfname))
628               {
629                 error("dicom_hdrinfo_extract: no such file or filename", warning);
630                 return (false);
631               }
632             
633             /*if ((store == NULL) || (store->to == NULL))
634               continue;*/
635             if ((fp = fopen_2(*pfname, "rb")) == NULL)
636               {       
637                 error("file pointer equals NULL", warning);
638                 return(false);
639               }
640             if ((type = dicom_preread_tests(fp)) == RUBBISH)
641               {
642                 format("error: file failed pre checks\n");
643                 return (false);
644               }
645             
646             fread(&group, sizeof(short), 1, fp);
647             short_swap((char *) &group);
648             
649             while (fp && !success && group != DCM_GROUPPIXEL)
650               {
651                 fread(&element, sizeof(short), 1, fp);
652                 short_swap((char *) &element);
653                 fread(&dblock, sizeof(int), 1, fp);
654                 word_swap((char *) &dblock);
655                 dblock_vr_conv(&dblock, fp);
656                 
657                 switch (DCM_MAKETAG(group, element))
658                   {
659                   case DCM_IDSERIESTIME:
660                     fread(temp, sizeof(char), dblock, fp);
661                     temp[dblock] = '\0';
662                     ser_time = get_dicom_time(temp);
663                     break;
664                     
665                   case DCM_IDIMAGETIME:
666                     fread(temp, sizeof(char), dblock, fp);
667                     temp[dblock] = '\0';
668                     im_time = get_dicom_time(temp);
669                     /*printf("image time    :%lf\n", im_time);*/
670                     success = true;
671                     break;
672                     
673                   case DCM_DLMITEM:
674                   case DCM_DLMITEMDELIMITATIONITEM:
675                   case DCM_DLMSEQUENCEDELIMITATIONITEM:
676                     break;
677                     
678                     
679                   default:
680                     for (k = 0; k < dblock; k++)
681                       fread(&junk, sizeof(char), 1, fp);
682                     break;
683                   }
684                 fread(&group, sizeof(short), 1, fp);
685                 short_swap((char *) &group);
686               }
687             
688             if (success)
689               {
690                 time = im_time - ser_time;
691                 times[j] =time;
692                 status = true;
693                 success = false;
694               }
695                     
696             fclose_2(fp, filename);
697             j++;
698           }
699         if (status)
700           {
701             /* seq->props = proplist_rm(seq->props, DYNSTIME);
702                seq->props = proplist_addifnp(seq->props, (void *) times, DYNSTIME, dynstimes_free, false);*/
703             seq->dim[3] = times[1];
704           }
705         
706         return (status);
707 }
708 
709 
710 Bool dicom_hdr_heartrate_extract(Sequence *seq, int end)
711 {
712   float *times=NULL, btime, phases, heartrate;
713   float beattime;
714   int type;
715   int psuccess = 0;
716   char temp[64];
717   unsigned char junk;
718   unsigned short group, element;
719   unsigned int dblock, tag1, tag2;
720   
721   Bool status = false, success = false;
722   FILE *fp = NULL;
723   char *pfname[1];
724   char filename[512];
725   char filename_temp[512];
726   int i, tend, k;
727   int j=0;
728 
729   strcpy(filename, "");
730   
731   /* mjs just temporarily, in futture extend with proper checks and stuff for 
732      all parameters */
733   tend = (end - seq->offset)+1;
734  
735   times = fvector_alloc(0, tend);
736  
737   for (i = (seq->offset); i <= end; i = i+seq->stride)
738     {
739       status = false;
740       parse_fname(seq->filename, filename_temp, i);
741       sprintf(filename, "%s%s", filename_temp, "");
742       *pfname = filename;
743       if (strchr(filename, '*') && !fname_resolve(NULL, filename, pfname))
744         {
745           error("dicom_hdrinfo_extract: no such file or filename", warning);
746           return (false);
747         }
748       
749       /*if ((store == NULL) || (store->to == NULL))
750         continue;*/
751       if ((fp = fopen_2(*pfname, "rb")) == NULL)
752         {       
753           error("file pointer equals NULL", warning);
754           return(false);
755         }
756       if ((type = dicom_preread_tests(fp)) == RUBBISH)
757         {
758           format("error: file failed pre checks\n");
759           return (false);
760         }
761       
762       fread(&group, sizeof(short), 1, fp);
763       short_swap((char *) &group);
764       
765       tag1 = 0x200117; /* ie, private group 2001,xx17 */
766       
767       while (fp && !success && group != DCM_GROUPPIXEL)
768         {
769           fread(&element, sizeof(short), 1, fp);
770           short_swap((char *) &element);
771           fread(&dblock, sizeof(int), 1, fp);
772           word_swap((char *) &dblock);
773           dblock_vr_conv(&dblock, fp);
774           
775           if (group == 0x2001)
776             {
777               format("0x2001 found");
778             }
779           tag2 = (group << 8) | (element & 0xFF); 
780           
781           if (tag2 == tag1) /* looks like it's 1017 anyway*/
782             {
783               fread(temp, sizeof(char), dblock, fp);
784               temp[dblock] = '\0';
785               sscanf(temp, "%f", &phases);
786               if (psuccess == 2)
787                 success = true;
788               else
789                 psuccess++;
790             }
791           else
792             {
793               
794               switch (DCM_MAKETAG(group, element))
795                 {/*
796                   case DCM_MAKETAG(0x0019, 0x1022):
797                   fread(temp, sizeof(char), dblock, fp);
798                   temp[dblock] = '\0';
799                   sscanf(temp, "%f", &btime);
800                   if (psuccess == 2)
801                     success = true;
802                   else
803                     psuccess++;
804                     break;*/
805                   /* mjs change no of phases tag - need to check*/
806                                 /* if 0x10 is right */
807                   
808                                 /* mjs change heart rate tag*/
809                 case DCM_MAKETAG(0x0018, 0x1088):
810                   fread(temp, sizeof(char), dblock, fp);
811                   temp[dblock] = '\0';
812                   sscanf(temp, "%f", &heartrate);
813                   /* if (psuccess == 1)*/
814                   success = true;
815                   /*else
816                     psuccess++;*/
817                   break;
818                   
819                 case DCM_DLMITEM:
820                 case DCM_DLMITEMDELIMITATIONITEM:
821                 case DCM_DLMSEQUENCEDELIMITATIONITEM:
822                   break;
823                   
824                   
825                 default:
826                   for (k = 0; k < dblock; k++)
827                     fread(&junk, sizeof(char), 1, fp);
828                   break;
829                 }
830             }
831           fread(&group, sizeof(short), 1, fp);
832           short_swap((char *) &group);
833         }
834       phases = 15;
835       if (phases == 1)
836         {
837           success = false;
838         }
839       if (success)
840         {
841           btime = 0.0;
842           beattime = 1.0 / (heartrate / 60);
843           times[j] = btime + (float) j *(beattime / (phases - 1));
844           success = false;
845           status = true;
846         }
847          
848       fclose_2(fp, filename);
849       j++;
850     }
851   if (status)
852     {
853       /*    seq->props = proplist_rm(seq->props, DYNSTIME);
854             seq->props = proplist_addifnp(seq->props, (void *) times, DYNSTIME, dynstimes_free, false);  */
855       seq->dim[3] = times[1];
856     }
857   
858   return (status);
859   
860 }
861 
862 
863 Bool dicom_hdr_patientdetails_extract()
864 {/* again, can assume that the details are the same for all images */
865   Bool success = false;
866   int type;
867   char *details = NULL;
868   unsigned char junk;
869   unsigned short group, element;
870   unsigned int dblock;
871   FILE *fp = NULL;
872   char *pfname[1];
873   char filename[512];
874   char filename_temp[512];
875   int  k;
876   Sequence *seq = (Sequence *)seq_get_current();
877 
878   strcpy(filename, "");
879    
880   parse_fname(seq->filename, filename_temp, seq->offset);
881   sprintf(filename, "%s%s", filename_temp, "");
882   *pfname = filename;
883   if (strchr(filename, '*') && !fname_resolve(NULL, filename, pfname))
884     {
885       error("dicom_hdrinfo_extract: no such file or filename", warning);
886       return (false);
887     }
888   
889   if ((fp = fopen_2(*pfname, "rb")) == NULL)
890     {
891       error("file pointer equals NULL", warning);
892       return(false);
893     }
894   if ((type = dicom_preread_tests(fp)) == RUBBISH)
895     {
896       format("error: file failed pre checks\n");
897       return (false);
898     }
899   
900   fread(&group, sizeof(short), 1, fp);
901   short_swap((char *) &group);
902   
903   details = (char *) ralloc(256*sizeof(char));
904   while (fp && !success && group != DCM_GROUPPIXEL)
905     {
906       fread(&element, sizeof(short), 1, fp);
907       short_swap((char *) &element);
908       fread(&dblock, sizeof(int), 1, fp);
909       word_swap((char *) &dblock);
910       dblock_vr_conv(&dblock, fp);
911       
912       switch (DCM_MAKETAG(group, element))
913         {
914         case DCM_MAKETAG(0x0008, 0x0020):
915           fread(details, sizeof(char), dblock, fp);
916           details[dblock] = '\0';
917           success = true;
918           break;
919           
920         case DCM_DLMITEM:
921         case DCM_DLMITEMDELIMITATIONITEM:
922         case DCM_DLMSEQUENCEDELIMITATIONITEM:
923           break;
924           
925           
926         default:
927           for (k = 0; k < dblock; k++)
928             fread(&junk, sizeof(char), 1, fp);
929           break;
930         }
931       fread(&group, sizeof(short), 1, fp);
932       short_swap((char *) &group);
933     
934     }
935   if (success != true)
936     cvector_free(details, 0);
937   else
938     {
939       seq->props = proplist_rm(seq->props, PAT_DATA);
940       seq->props = proplist_addifnp(seq->props, (void *) details, PAT_DATA, rfree, false);
941     }
942   
943   fclose_2(fp, filename);
944   return (success);
945   
946 }
947 
948 
949 Bool dicom_hdr_voxelscale_extract()
950 {/*logically only need to look at the first image. */
951   float xsize, ysize, zsize;
952   int type;
953   int psuccess = 0;
954   char dimension[64];
955   unsigned char junk;
956   unsigned short group, element;
957   unsigned int dblock;
958   Bool success = false;
959   FILE *fp = NULL;
960   char *pfname[1];
961   char filename[512];
962   char filename_temp[512];
963   int k;
964   Sequence * seq = (Sequence *)seq_get_current();
965 
966 
967   strcpy(filename, "");
968   
969   parse_fname(seq->filename, filename_temp, seq->offset);
970   sprintf(filename, "%s%s", filename_temp, "");
971   *pfname = filename;
972   if (strchr(filename, '*') && !fname_resolve(NULL, filename, pfname))
973     {
974       error("dicom_hdrinfo_extract: no such file or filename", warning);
975       return (false);
976     }
977   
978   if ((fp = fopen_2(*pfname, "rb")) == NULL)
979     {  
980       error("file pointer equals NULL", warning);
981       return(false);
982     }
983   xsize = ysize = zsize = 1.0;
984   if ((type = dicom_preread_tests(fp)) == RUBBISH)
985     {
986       format("error: file failed pre checks\n");
987       return (false);
988     }
989   
990   fread(&group, sizeof(short), 1, fp);
991   short_swap((char *) &group);
992   
993   while (fp && group != DCM_GROUPPIXEL && psuccess < 3)
994     {
995       
996       fread(&element, sizeof(short), 1, fp);
997       short_swap((char *) &element);
998       fread(&dblock, sizeof(int), 1, fp);
999       word_swap((char *) &dblock);
1000       dblock_vr_conv(&dblock, fp);
1001       
1002       switch (DCM_MAKETAG(group, element))
1003         {
1004         case DCM_ACQSLICESPACING:
1005           fread(dimension, sizeof(char), dblock, fp);
1006           dimension[dblock] = '\0';
1007           sscanf(dimension, "%f", &zsize);
1008           psuccess++;
1009           break;
1010           
1011         case DCM_IMGPIXELSPACING:
1012           fread(dimension, sizeof(char), dblock, fp);
1013           dimension[dblock] = '\0';
1014           sscanf(dimension, "%f\\%f", &xsize, &ysize);
1015           psuccess++;
1016           break;
1017           
1018         case DCM_DLMITEM:
1019         case DCM_DLMITEMDELIMITATIONITEM:
1020         case DCM_DLMSEQUENCEDELIMITATIONITEM:
1021           break;
1022           
1023           
1024         default:
1025           for (k = 0; k < dblock; k++)
1026             fread(&junk, sizeof(char), 1, fp);
1027           break;
1028         }
1029       fread(&group, sizeof(short), 1, fp);
1030       short_swap((char *) &group);
1031     }
1032   
1033   if (psuccess == 2)
1034     {
1035       /*     
1036              iscale = vec3_alloc();
1037              iscale->el[0] = xsize;
1038              iscale->el[1] = ysize;
1039              iscale->el[2] = zsize;
1040       */
1041 
1042       seq->dim[0] = xsize;
1043       seq->dim[1] = ysize;
1044       seq->dim[2] = zsize;
1045 
1046       success = true;
1047       /*
1048         seq->props = proplist_rm(seq->props, VOXELS);
1049         seq->props = proplist_addifnp(seq->props, iscale, VOXELS, vec3_free, false);
1050       */
1051     }
1052   
1053   
1054   fclose_2(fp, filename);
1055   
1056   return success;
1057 }
1058 
1059 static Imrect *dicom_read_multiformat_image(const char *pathname, FILE * fp_img, int scale_factor_flag)
1060 {
1061         int get_swapping_ts();
1062         Imrect *imrect = NULL, *imrect2 = NULL;
1063         Imregion imregion;
1064         Vartype new_vtype;
1065         float range;
1066         float rinter, rslope;
1067         float sinter, sslope; /*mjs added 25/11/05 to give functionality for loading scale slope etc */
1068         float pix;
1069         int i, j, k;
1070         unsigned char junk;
1071         char dimension[64];
1072         unsigned short abits;
1073         unsigned short sbits;
1074         unsigned short sign;
1075         unsigned short rows;
1076         unsigned short cols;
1077         unsigned short group;
1078         unsigned short element;
1079         unsigned int dblock;
1080         float xsize, ysize, zsize;
1081         Sequence *seq = (Sequence *) seq_get_current();
1082        
1083         rinter = 0.0;
1084         rslope = 1.0;
1085         xsize = ysize = zsize = 1.0;
1086 
1087         fread(&group, sizeof(short), 1, fp_img);
1088         short_swap((char *) &group);
1089         fread(&element, sizeof(short), 1, fp_img);
1090         short_swap((char *) &element);
1091 
1092         while (fp_img && DCM_MAKETAG(group, element) != DCM_PXLPIXELDATA)
1093         {
1094                 fread(&dblock, sizeof(int), 1, fp_img);
1095                 word_swap((char *) &dblock);
1096                 dblock_vr_conv(&dblock, fp_img);
1097 
1098                 if (DCM_MAKETAG(group, element) == DCM_PXLPIXELDATA)
1099                         break;
1100 
1101                 switch (DCM_MAKETAG(group, element))
1102                 {
1103                         case DCM_ACQSLICESPACING:
1104                                 fread(dimension, sizeof(char), dblock, fp_img);
1105                                 dimension[dblock] = '\0';
1106                                 sscanf(dimension, "%f", &zsize);
1107                                 break;
1108 
1109                         case DCM_IMGPIXELSPACING:
1110                                 fread(dimension, sizeof(char), dblock, fp_img);
1111                                 dimension[dblock] = '\0';
1112                                 sscanf(dimension, "%f\\%f", &xsize, &ysize);
1113                                 break;
1114 
1115                         case DCM_IMGRESCALEINTERCEPT:
1116                                 fread(dimension, sizeof(char), dblock, fp_img);
1117                                 dimension[dblock] = '\0';
1118                                 sscanf(dimension, "%f", &rinter);
1119                                 break;
1120 
1121                                 /*case DCM_IMGRESCALESLOPE:*/
1122                 case DCM_MAKETAG(0x0028,0x1053):        
1123                                 fread(dimension, sizeof(char), dblock, fp_img);
1124                                 dimension[dblock] = '\0';
1125                                 sscanf(dimension, "%f", &rslope);
1126                                 break;
1127                                 /*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.*/
1128                         case DCM_MAKETAG(0x2005,0x100d):
1129                           fread(&sinter, sizeof(float), 1, fp_img);
1130                           break;
1131                         case DCM_MAKETAG(0x2005,0x100e):
1132                     
1133                           fread(&sslope, sizeof(float), 1, fp_img);
1134                           break;
1135                                 
1136                         case DCM_IMGROWS:
1137                                 fread(&rows, sizeof(short), 1, fp_img);
1138                                 short_swap((char *) &rows);
1139                                 break;
1140 
1141                         case DCM_IMGCOLUMNS:
1142                                 fread(&cols, sizeof(short), 1, fp_img);
1143                                 short_swap((char *) &cols);
1144                                 break;
1145 
1146                         case DCM_IMGPIXELREPRESENTATION:
1147                                 fread(&sign, sizeof(short), 1, fp_img);
1148                                 short_swap((char *) &sign);
1149                                 break;
1150 
1151                         case DCM_IMGBITSSTORED:
1152                                 fread(&sbits, sizeof(short), 1, fp_img);
1153                                 short_swap((char *) &sbits);
1154                                 break;
1155 
1156                         case DCM_IMGBITSALLOCATED:
1157                                 fread(&abits, sizeof(short), 1, fp_img);
1158                                 short_swap((char *) &abits);
1159                                 break;
1160 
1161                         case DCM_DLMITEM:
1162                         case DCM_DLMITEMDELIMITATIONITEM:
1163                         case DCM_DLMSEQUENCEDELIMITATIONITEM:
1164                                 break;
1165 
1166                         default:
1167                                 for (i = 0; i < dblock; i++)
1168                                         fread(&junk, sizeof(char), 1, fp_img);
1169                                 break;
1170                 }
1171                 fread(&group, sizeof(short), 1, fp_img);
1172                 short_swap((char *) &group);
1173                 fread(&element, sizeof(short), 1, fp_img);
1174                 short_swap((char *) &element);
1175         }
1176 
1177         if (fp_img)
1178         {
1179                 /* 
1180                  * move to start of image data depending on VR mode
1181                  */
1182                 fread(&dblock, sizeof(int), 1, fp_img);
1183                 word_swap((char *) &dblock);
1184                 dblock_vr_conv(&dblock, fp_img);
1185 
1186                 imregion.lx = 0;
1187                 imregion.ux = cols;
1188                 imregion.ly = 0;
1189                 imregion.uy = rows;
1190 
1191                 if (sign == 0)
1192                 {
1193                         if (abits == 16 || abits == 12)
1194                                 new_vtype = ushort_v;
1195                         else
1196                                 new_vtype = uchar_v;
1197                 } else
1198                 {
1199                         if (abits == 16 || abits == 12)
1200                                 new_vtype = short_v;
1201                         else
1202                                 new_vtype = char_v;
1203                 }
1204                 if (abits == 12)
1205                         imregion.ux = 3.0 * imregion.ux / 4;
1206                 cols = imregion.ux;
1207 
1208                 imrect = im_alloc(rows, cols, &imregion, (Vartype) new_vtype);
1209 
1210                 if (!fread_imrect_data(imrect, fp_img, pathname))
1211                 {
1212                         im_free(imrect);
1213                         imrect = NULL;
1214                         return (NULL);
1215                 }
1216                 if (abits == 12)
1217                 {
1218                         imrect2 = im_dicom_conv(imrect);
1219                         im_free(imrect);
1220                         imrect = imrect2;
1221                 }
1222                 if (imrect != NULL && get_swapping_ts())
1223                 {
1224                         im_endian_inplace(imrect);
1225                 }
1226 
1227                 range = pow(2.0, sbits) - 1;
1228 
1229 
1230                 imrect2 = im_cast(imrect, float_v);
1231                 im_free(imrect);
1232                 imrect = imrect2;
1233                 if (rslope != 0.0 || rinter != 0.0)
1234                 {
1235                         for (j = imrect->region->ly; j < imrect->region->uy; j++)
1236                                 for (k = imrect->region->lx; k < imrect->region->ux; k++)
1237                                 {
1238                                         pix = im_get_pixf(imrect, j, k);
1239                                         if (scale_factor_flag == 0)
1240                                           {
1241                                             pix = pix * rslope + rinter;
1242                                           }
1243                                         if (scale_factor_flag == 1)
1244                                           {
1245                                             pix = (pix-sinter)/sslope;
1246                                           }
1247                                         im_put_pixf(pix, imrect, j, k);
1248                                 }
1249                 }
1250         
1251                   if (seq != NULL)
1252                   {
1253                           seq->dim[0] = xsize;
1254                           seq->dim[1] = ysize;
1255                           seq->dim[2] = zsize;
1256                   }
1257                 
1258         }
1259 
1260         return imrect;
1261 }
1262 
1263 
1264 
1265 Imrect *dicom_read_image(char *pathname, int scale_factor_flag)
1266 {
1267         FILE *fp_img;
1268         Imrect *im;
1269 
1270         if ((fp_img = fopen_2(pathname, "rb")) == NULL)
1271                 return (NULL);
1272 
1273         switch (dicom_preread_tests(fp_img))
1274         {
1275                 case PART10:
1276                         im = dicom_read_multiformat_image(pathname, fp_img, scale_factor_flag);
1277                         fclose_2(fp_img, pathname);
1278                         break;
1279                 case NONPART10:
1280                         im = dicom_read_multiformat_image(pathname, fp_img, scale_factor_flag);
1281                         fclose_2(fp_img, pathname);
1282                         break;
1283                 default:
1284                         format("error: file failed pre checks\n");
1285                         fclose_2(fp_img, pathname);
1286                         im = NULL;
1287                         break;
1288         }
1289 
1290         return (im);
1291 }
1292 
1293 
1294 /*
1295   recover information from dicom headers ref'd in seq using hdrextract_func
1296   which should store info in image propslist.  hdrextract_func should be
1297   passed fileptr and imrect and should return Bool success. 
1298 */
1299 /*
1300 Bool dicom_hdrinfo_extract(Sequence * seq, void (*hdrextract_func))
1301 {
1302         List *store = (List *) get_seq_start_el(seq);
1303         Bool status = false, success = true;
1304         FILE *fp = NULL;
1305         char *pfname[1];
1306         char filename[512];
1307         char temp[512];
1308         int i;
1309         int end = get_end_frame(seq);
1310 
1311         strcpy(filename, "");
1312 
1313         for (i = (seq->offset); i <= end; i = i+seq->stride)
1314         {
1315                 parse_fname(seq->filename, temp, i);
1316                 sprintf(filename, "%s%s", temp, "");
1317                 *pfname = filename;
1318                 if (strchr(filename, '*') && !fname_resolve(NULL, filename, pfname))
1319                 {
1320                         error("dicom_hdrinfo_extract: no such file or filename", warning);
1321                         return (false);
1322                 }
1323 
1324                 if ((store == NULL) || (store->to == NULL))
1325                         continue;
1326                 if ((fp = fopen_2(*pfname, "rb")) == NULL)
1327                         continue;
1328                 status = ((Bool(*)())hdrextract_func) (fp, (Imrect *) (store->to));
1329                 fclose_2(fp, filename);
1330                 if (!status)
1331                         success = status;
1332                 store = store->next;
1333         }
1334 
1335         return (success);
1336 }
1337 */
1338 

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