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

Linux Cross Reference
Tina4/src/file/nema/nema_io.c

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

  1 /**@(#)Functions to read/write ACR-NEMA image files
  2 */
  3 
  4 #include <stdio.h>
  5 #include <string.h>
  6 #include <tina/file.h>
  7 #include <tina/sys.h>
  8 #include <tina/sys_types.h>
  9 #include <tina/sysfuncs.h>
 10 #include <tina/math.h>
 11 #include <tina/mathfuncs.h>
 12 #include <tina/temporal.h>
 13 
 14 #define VOXELS    450
 15 #define DYNSTIME  451
 16 #define TE_DATA   453
 17 #define PAT_DATA  454
 18 
 19 
 20 
 21 extern Bool  fclose_2(FILE * stream, const char *pathname);
 22 extern FILE *fopen_2(const char *pathname, const char *mode);
 23 extern Bool  fread_imrect_data(const Imrect * imrect, FILE * stream, const char *pathname);
 24 extern Bool  fwrite_imrect_data(const Imrect * imrect, FILE * stream, const char *pathname);
 25 
 26 void  dynstimes_free(float *times)
 27 {
 28   fvector_free(times, 0);
 29 }
 30 
 31 
 32 /* Create an imrect and read an image file called 'pathname' into it.
 33  * Open file, read data and close file. On failure, give error message
 34  * and return NULL. */
 35 
 36 static int valid_group(short group)
 37 {
 38         if(group == 0x0000 || group == 0x0008 || group == 0x0010
 39         || group == 0x0020 || group == 0x0018 || group == 0x0020
 40         || group == 0x0028 || group == 0x4000 || group == 0x4000
 41         || group == 0x7FE0 )
 42            return(1);
 43         else
 44            return(0);
 45 } 
 46 
 47 Imrect *im_nema_conv(Imrect *im)
 48 {
 49     Imrect *im2;
 50     Imregion roi;
 51     unsigned char *row1;
 52     unsigned char *row2;
 53     int i,j,k;
 54     int lx, ux, ly, uy;
 55 
 56     if (im == NULL)
 57       return(NULL);
 58 
 59     roi = *(im->region);
 60     roi.ux = roi.ux*4.0/3;
 61     im2 = im_alloc(im->height,roi.ux,&roi,short_v);
 62     lx = roi.lx;
 63     ux = roi.ux;
 64     ly = roi.ly;
 65     uy = roi.uy;
 66 
 67     for (i = ly; i < uy; ++i)
 68     {
 69         row1 = IM_ROW(im,i);
 70         row2 = IM_ROW(im2,i);
 71         for (k = 2*lx, j = 2*lx; k < 2*ux; j+=3)
 72         {
 73             row2[k++] = row1[j];
 74             row2[k++] = row1[j+1]%16;
 75             row2[k++] = row1[j+1]/16 + (row1[j+2]%16)*16;
 76             row2[k++] = row1[j+2]/16;
 77         }
 78     }
 79 
 80     return(im2);
 81 }
 82 
 83 void im_endian_inplace(Imrect *imrect)
 84 {
 85     unsigned int size = var_size(imrect->vtype);
 86     Imregion roi;
 87     int y;
 88     int lx, ux, ly, uy;
 89 
 90     roi = *(imrect->region);
 91     lx = roi.lx;
 92     ux = roi.ux;
 93     ly = roi.ly;
 94     uy = roi.uy;
 95 
 96 
 97     for(y = ly; y < uy; y++)
 98     {
 99         void*  row = NULL;
100         char   temp[40];
101         int    x;
102         unsigned int i;
103 
104         row = IM_ROW(imrect,y);
105 
106         for(x = lx; x < ux; x++)
107         {
108             for(i = 0; i < size; i++)
109                 temp[i] = *(((char *)row) + (x*size+i));
110             for(i = 1; i <= size; i++)
111                 *(((char *)row) + (x*size + size-i)) = temp[i-1];
112         }
113     }
114 }
115 
116 
117 Bool   nema_hdr_TE_extract(FILE *fp, Imrect *im)
118 {
119   Bool                swapit = false;
120   Bool                success = false;
121   float              *TE;
122   int                 i, j, k;
123   char               *buffer = NULL;
124   char                temp[64];
125   unsigned char       junk;
126   unsigned short      group, element;
127   unsigned short      dblock;
128   unsigned short      dblock1;
129 
130   set_swapping_ts(0);
131   TE = (float *)ralloc(sizeof(float));
132   while (fp && fread(&group, sizeof(short), 1, fp) && !success)
133     { 
134       if(!swapit)
135         {
136           if(!valid_group(group))
137             {
138               set_swapping_ts(1);
139               short_swap((char *)&group);
140               if (valid_group(group)) 
141                 swapit = true;
142               else    
143                 {
144                   short_swap((char *)&group);
145                   set_swapping_ts(0);
146                 }
147             }
148         }
149       else
150         short_swap((char *)&group);
151  
152       fread(&element, sizeof(short), 1, fp);
153       short_swap((char *)&element);
154       fread(&dblock,  sizeof(short), 1, fp);
155       short_swap((char *)&dblock);
156       fread(&dblock1, sizeof(short), 1, fp);
157       short_swap((char *)&dblock1);
158 
159       if (dblock1 != 0)
160         {
161           format("error: non-standard block size\n");
162           return(false);
163         }
164       else if (group == 0x0018 && element == 0x0081)
165         {
166           fread(temp, sizeof(char), dblock, fp);
167           temp[dblock] = '\0';
168           sscanf(temp, "%f", TE);
169           *TE = *TE/100.0;
170           success = true;
171         }
172       else 
173         for (i = 0; i < dblock; i++) 
174           fread(&junk, sizeof(char), 1, fp);
175 
176       if (success) break;
177     }
178 
179   if (success)
180     {
181       im->props = proplist_rm(im->props, TE_DATA);
182       im->props = proplist_addifnp(im->props, (void *)TE, TE_DATA, rfree, false);
183     }
184 
185   return (success);
186 }
187 
188 
189 Bool   nema_hdr_dynstimes_extract(FILE *fp, Imrect *im)
190 {
191   Bool                swapit = false;
192   Bool                success = false;
193   float              *times;
194   int                 i, j, k, stime;
195   int                 dynstudy;
196   char               *buffer = NULL;
197   char                temp[64], v;
198   unsigned char       junk;
199   unsigned short      group, element;
200   unsigned short      dblock;
201   unsigned short      dblock1;
202 
203   set_swapping_ts(0);
204   while (fp && fread(&group, sizeof(short), 1, fp) && !success)
205     { 
206       if(!swapit)
207         {
208           if(!valid_group(group))
209             {
210               set_swapping_ts(1);
211               short_swap((char *)&group);
212               if (valid_group(group)) 
213                 swapit = true;
214               else    
215                 {
216                   short_swap((char *)&group);
217                   set_swapping_ts(0);
218                 }
219             }
220         }
221       else
222         short_swap((char *)&group);
223  
224       fread(&element, sizeof(short), 1, fp);
225       short_swap((char *)&element);
226       fread(&dblock,  sizeof(short), 1, fp);
227       short_swap((char *)&dblock);
228       fread(&dblock1, sizeof(short), 1, fp);
229       short_swap((char *)&dblock1);
230 
231       if (dblock1 != 0)
232         {
233           format("error: non-standard block size\n");
234           return(false);
235         }
236       else if (group == 0x0019 && element == 0x1015)
237         {
238           fread(temp, sizeof(char), dblock, fp);
239           temp[dblock] = '\0';
240           sscanf(temp, "%d", &dynstudy);
241           if (!dynstudy)
242             {
243               return(false);
244             }
245         }
246       else if (group == 0x0019 && element == 0x101B)
247         {
248           fread(temp, sizeof(char), dblock, fp);
249           temp[dblock] = '\0';
250           sscanf(temp, "%d", &stime);
251         } 
252       else if (group == 0x0019 && element == 0x1021)
253         {
254           buffer = (char *)ralloc(dblock*sizeof(char));
255           times = fvector_alloc(0, stime);
256           fread(buffer, sizeof(char), dblock, fp);
257 
258           j = k = 0;
259           for (i = 0; i < stime; i++)
260             {
261               while (buffer[j] != '\\')
262                 temp[k++] = buffer[j++];
263               temp[k] = '\0';
264               sscanf(temp, "%f", &(times[i]));
265               k = 0; j++;
266             }
267           rfree(buffer);
268           success = true;
269           break;
270         } 
271       else 
272         for (i = 0; i < dblock; i++) 
273           fread(&junk, sizeof(char), 1, fp);
274     }
275 
276   if (success)
277     {
278       im->props = proplist_rm(im->props, DYNSTIME);
279       im->props = proplist_addifnp(im->props, (void *)times, DYNSTIME, dynstimes_free, false);
280     }
281 
282   return (success);
283 }
284 
285 
286 Bool   nema_hdr_heartrate_extract(FILE *fp, Imrect *im)
287 {
288   Bool                swapit = false;
289   Bool                success = false;
290   float              *times, btime, phases, heartrate;
291   float               beattime;
292   int                 i, j, k;
293   int                 psuccess = 0;
294   char               *buffer = NULL;
295   char                temp[64], v;
296   unsigned char       junk;
297   unsigned short      group, element;
298   unsigned short      dblock;
299   unsigned short      dblock1;
300 
301   set_swapping_ts(0);
302   while (fp && fread(&group, sizeof(short), 1, fp) && !success)
303     { 
304       if(!swapit)
305         {
306           if(!valid_group(group))
307             {
308               set_swapping_ts(1);
309               short_swap((char *)&group);
310               if (valid_group(group)) 
311                 swapit = true;
312               else    
313                 {
314                   short_swap((char *)&group);
315                   set_swapping_ts(0);
316                 }
317             }
318         }
319       else
320         short_swap((char *)&group);
321  
322       fread(&element, sizeof(short), 1, fp);
323       short_swap((char *)&element);
324       fread(&dblock,  sizeof(short), 1, fp);
325       short_swap((char *)&dblock);
326       fread(&dblock1, sizeof(short), 1, fp);
327       short_swap((char *)&dblock1);
328 
329       if (dblock1 != 0)
330         {
331           format("error: non-standard block size\n");
332           return(false);
333         }
334       else if (group == 0x0019 && element == 0x1022)
335         {
336           fread(temp, sizeof(char), dblock, fp);
337           temp[dblock] = '\0';
338           sscanf(temp, "%f", &btime);
339           if (psuccess == 2) success = true;
340           else psuccess++;
341         }
342       else if (group == 0x0019 && element == 0x1069)
343         {
344           fread(temp, sizeof(char), dblock, fp);
345           temp[dblock] = '\0';
346           sscanf(temp, "%f", &phases);
347           if (psuccess == 2) success = true;
348           else psuccess++;
349         }
350        else if (group == 0x0019 && element == 0x106A)
351         {
352           fread(temp, sizeof(char), dblock, fp);
353           temp[dblock] = '\0';
354           sscanf(temp, "%f", &heartrate);
355           if (psuccess == 2) success = true;
356           else psuccess++;
357         }
358       else 
359         for (i = 0; i < dblock; i++) 
360           fread(&junk, sizeof(char), 1, fp);
361     }
362 
363   if (phases == 1) success = false;
364   if (success)
365     {
366       times = fvector_alloc(0, phases);
367       beattime = 1.0/(heartrate/60);
368       for (i = 0; i < phases; i++)
369         times[i] = btime + (float)i*(beattime/(phases-1));
370       im->props = proplist_rm(im->props, DYNSTIME);
371       im->props = proplist_addifnp(im->props, (void *)times, DYNSTIME, rfree, false);
372     }
373 
374   return (success);
375 }
376 
377 
378 Bool  nema_hdr_patientdetails_extract(FILE *fp, Imrect *im)
379 {
380   Bool                swapit = false;
381   Bool                psuccess = false;
382   int                 i;
383   char               *details = NULL;
384   unsigned char       junk;
385   unsigned short      group, element;
386   unsigned short      dblock;
387   unsigned short      dblock1;
388 
389   details = (char *)cvector_alloc(0, 256);
390   set_swapping_ts(0);
391   while (fp && fread(&group, sizeof(short), 1, fp) && !psuccess)
392     { 
393       if(!swapit)
394         {
395           if(!valid_group(group))
396             {
397               set_swapping_ts(1);
398               short_swap((char *)&group);
399               if (valid_group(group)) 
400                 swapit = true;
401               else    
402                 {
403                   short_swap((char *)&group);
404                   set_swapping_ts(0);
405                 }
406             }
407         }
408       else
409         short_swap((char *)&group);
410  
411       fread(&element, sizeof(short), 1, fp);
412       short_swap((char *)&element);
413       fread(&dblock,  sizeof(short), 1, fp);
414       short_swap((char *)&dblock);
415       fread(&dblock1, sizeof(short), 1, fp);
416       short_swap((char *)&dblock1);
417 
418       if (dblock1 != 0)
419         {
420           format("error: non-standard block size\n");
421           cvector_free(details, 0);
422           return(false);
423         }
424       else if (group == 0x0008 && element == 0x0020)
425         {
426           fread(details, sizeof(char), dblock, fp);
427           details[dblock] = '\0';
428           psuccess = true;
429         }
430       else 
431         for (i = 0; i < dblock; i++) 
432           fread(&junk, sizeof(char), 1, fp);
433     }
434 
435   if (psuccess != true)
436                 cvector_free(details, 0);
437         else
438                 {
439       im->props = proplist_rm(im->props, PAT_DATA);
440       im->props = proplist_addifnp(im->props, (void *)details, PAT_DATA, 
441                                                                                                                                          rfree, false);
442                 }
443 
444   return (psuccess);
445 }
446 
447 
448 Bool  nema_hdr_voxelscale_extract(FILE *fp, Imrect *im)
449 {
450   Bool                swapit = false;
451         Vec3               *iscale;
452   float               xsize,ysize,zsize;
453   float               zgap;
454   int                 psuccess = 0;
455   int                 i;
456   char               *details = NULL;
457   char                dimension[64];
458   unsigned char       junk;
459   unsigned short      group, element;
460   unsigned short      dblock;
461   unsigned short      dblock1;
462 
463   details = (char *)cvector_alloc(0, 256);
464   set_swapping_ts(0);
465   while (fp && fread(&group, sizeof(short), 1, fp) && psuccess < 3)
466     { 
467       if(!swapit)
468         {
469           if(!valid_group(group))
470             {
471               set_swapping_ts(1);
472               short_swap((char *)&group);
473               if (valid_group(group)) 
474                 swapit = true;
475               else    
476                 {
477                   short_swap((char *)&group);
478                   set_swapping_ts(0);
479                 }
480             }
481         }
482       else
483         short_swap((char *)&group);
484  
485       fread(&element, sizeof(short), 1, fp);
486       short_swap((char *)&element);
487       fread(&dblock,  sizeof(short), 1, fp);
488       short_swap((char *)&dblock);
489       fread(&dblock1, sizeof(short), 1, fp);
490       short_swap((char *)&dblock1);
491 
492       if (dblock1 != 0)
493         {
494           format("error: non-standard block size\n");
495           return(false);
496         }
497       else if (group == 0x0018 && element == 0x0050)
498         {
499           fread(dimension, sizeof(char),dblock,fp);
500           dimension[dblock] = '\0';
501           sscanf(dimension,"%f",&zsize);
502                                         psuccess++;
503         }
504       else if (group == 0x0021 && element == 0x1221)
505         {
506           fread(dimension, sizeof(char),dblock,fp);
507           dimension[dblock] = '\0';
508           sscanf(dimension,"%f",&zgap);
509                                         psuccess++;
510         }
511                         else if(group == 0x0028 && element == 0x0030)
512                                 {
513                                         fread(dimension, sizeof(char),dblock,fp);
514                                         dimension[dblock] = '\0';
515                                         sscanf(dimension,"%f\\%f",&xsize,&ysize);
516                                         psuccess++;
517                                 }
518  
519       else 
520         for (i = 0; i < dblock; i++) 
521           fread(&junk, sizeof(char), 1, fp);
522     }
523 
524   if (psuccess == 3)
525                 {
526       iscale = vec3_alloc();
527       iscale->el[0] = xsize;
528       iscale->el[1] = ysize;
529       iscale->el[2] = zsize + zgap;
530 
531       im->props = proplist_rm(im->props, VOXELS);
532       im->props = proplist_addifnp(im->props, iscale,
533                                                                                                                                          VOXELS, vec3_free, false);
534                 }
535 
536   return ((Bool)(psuccess == 3));
537 }
538 
539 
540 Imrect *nema_read_image(const char *pathname)
541 {
542   Imrect *imrect = NULL, *imrect2 = NULL, *fim = NULL;
543   Imregion imregion;
544   FILE                    *fp_img=fopen_2(pathname, "rb");
545   Vartype                new_vtype;
546   float        fpmin, fpmax, range;
547   float                   val, pix;
548   int                      i, j, k;
549   short                   endian=0;
550   unsigned char               junk;
551   char               dimension[64];
552   unsigned short              abits;
553   unsigned short              sbits;
554   unsigned short              sign;
555   unsigned short              rows;
556   unsigned short              cols;
557   unsigned short             group;
558   unsigned short           element;
559   unsigned short            dblock;
560   unsigned short           dblock1;
561   float          xsize,ysize,zsize;
562   float                       zgap;
563   Vec3                     *iscale;
564 
565   set_swapping_ts(0);
566   while (fp_img&&fread(&group, sizeof(short),1,fp_img))
567     {
568       if(!endian)
569         {
570           if(valid_group(group)!=1)
571             {
572               set_swapping_ts(1);
573               short_swap((char *)&group);
574               if (valid_group(group) == 1) endian = 1;
575               else    
576                 {
577                   short_swap((char *)&group);
578                   set_swapping_ts(0);
579                 }
580             }
581         }
582       else
583         {
584           short_swap((char *)&group);
585         } 
586       fread(&element, sizeof(short),1,fp_img);
587       short_swap((char *)&element);
588       fread(&dblock,  sizeof(short),1,fp_img);
589       short_swap((char *)&dblock);
590       fread(&dblock1, sizeof(short),1,fp_img);
591       short_swap((char *)&dblock1);
592 
593       if (group == 0x7FE0 && element == 0x0010)
594         break; /* found the image data */
595       if (dblock1!=0)
596         {
597           format("error: non-standard block size\n");
598           fclose_2(fp_img,pathname); 
599           return(NULL);
600         }
601       else if (group == 0x0018 && element == 0x0050)
602         {
603           fread(dimension, sizeof(char),dblock,fp_img);
604           dimension[dblock] = '\0';
605           sscanf(dimension,"%f",&zsize);
606         }
607       else if (group == 0x0021 && element == 0x1221)
608         {
609           fread(dimension, sizeof(char),dblock,fp_img);
610           dimension[dblock] = '\0';
611           sscanf(dimension,"%f",&zgap);
612         }
613       else if (group == 0x0029 && element == 0x1010)
614         {
615           fread(dimension, sizeof(char), dblock, fp_img);
616           dimension[dblock] = '\0';
617           sscanf(dimension, "%f", &fpmin);
618         }
619       else if (group == 0x0029 && element == 0x1011)
620         {
621           fread(dimension, sizeof(char), dblock, fp_img);
622           dimension[dblock] = '\0';
623           sscanf(dimension, "%f", &fpmax);
624         }
625 
626       else if (group == 0x0028)
627         {
628           if(element == 0x0010)
629             {
630               fread(&rows, sizeof(short),1,fp_img);
631               short_swap((char *)&rows);
632             }
633           else if(element == 0x0011)
634             {
635               fread(&cols, sizeof(short),1,fp_img);
636               short_swap((char *)&cols);
637             }
638           else if(element == 0x0103)
639             {
640               fread(&sign, sizeof(short),1,fp_img);
641               short_swap((char *)&sign);
642             }
643           else if(element == 0x0101)
644             {
645               fread(&sbits, sizeof(short),1,fp_img);
646               short_swap((char *)&sbits);
647             }
648           else if(element == 0x0100)
649             {
650               fread(&abits, sizeof(short),1,fp_img);
651               short_swap((char *)&abits);
652             }
653           else if(element == 0x0030)
654             {
655               fread(dimension, sizeof(char),dblock,fp_img);
656               dimension[dblock] = '\0';
657               sscanf(dimension,"%f\\%f",&xsize,&ysize);
658             }
659           else for (i=0;i<dblock;i++) 
660             fread(&junk, sizeof(char),1,fp_img);
661         }   
662       else 
663         for (i=0;i<dblock;i++) fread(&junk, sizeof(char),1,fp_img);
664     }
665 
666     
667   if (fp_img)
668     {
669       imregion.lx = 0;
670       imregion.ux = cols;
671       imregion.ly = 0;
672       imregion.uy = rows;
673         
674       format("x_dim=%i ",cols);
675       format("y_dim=%i\n",rows);
676       if (sign == 0)
677         {
678           if (abits == 16 || abits == 12) 
679             {
680               new_vtype = ushort_v;
681             }
682           else new_vtype = uchar_v;
683         }
684       else
685         { 
686           if (abits == 16 || abits == 12) 
687             {
688               new_vtype = short_v;
689             }
690           else new_vtype = char_v;
691         }
692       if(abits == 12) imregion.ux = 3.0*imregion.ux/4;
693       cols = imregion.ux;
694 
695       imrect = im_alloc(rows, cols, &imregion, (Vartype)new_vtype);
696 
697       if (!fread_imrect_data(imrect, fp_img, pathname))
698         {
699           im_free(imrect);
700           imrect = NULL;
701           fclose_2(fp_img, pathname);
702           return(NULL);
703         }
704       fclose_2(fp_img, pathname);
705       if (abits == 12)
706         {
707           imrect2 = im_nema_conv(imrect);
708           im_free(imrect);
709           imrect = imrect2;
710         }
711       if (imrect!=NULL&&endian)
712         {
713           im_endian_inplace(imrect);
714         }
715 
716       format("scale min %.2f  scale max %.2f\n", fpmin, fpmax);
717       range = pow(2.0, sbits) - 1;
718       fim = im_alloc(imrect->height, imrect->width, imrect->region, 
719                      (Vartype)float_v);
720       for(j = imrect->region->ly; j < imrect->region->uy; j++)
721         for(k = imrect->region->lx; k < imrect->region->ux; k++)
722           {
723             if ((pix = im_get_pixf(imrect, j, k)) != 0.0)
724                val = ((fpmax-fpmin)*pix)/(float)range + fpmin;
725             else
726                val = 0.0; 
727             im_put_pixf(val, fim, j, k);    
728           }   
729 
730       iscale = vec3_alloc();
731       iscale->el[0] = xsize;
732       iscale->el[1] = ysize;
733       iscale->el[2] = zsize + zgap;
734       fim->props = proplist_addifnp(fim->props, iscale,
735                                        VOXELS, vec3_free, false);
736 
737     }
738   
739   im_free(imrect);
740   return fim;
741 }
742 
743 
744 /*
745   recover information from nema headers ref'd in seq using hdrextract_func
746   which should store info in image propslist.  hdrextract_func should be
747   passed fileptr and imrect and should return Bool success. 
748 */
749 Bool    nema_hdrinfo_extract(Sequence *seq, void (*hdrextract_func))
750 {
751   List  *store = (List *)get_start_el();
752   Bool     status = false, success = true ;
753   FILE    *fp = NULL;
754   char    *pfname[1];
755   char     filename[512];
756   char     temp[512];
757   int      i;
758 
759   strcpy(filename, "");
760   strcpy(temp, "");
761  
762   for(i = (seq->seq_start); i <= seq->seq_end; i++)
763     {
764       parse_fname(seq->filename, temp, i);
765       sprintf(filename, "%s%s", temp, ".ani");
766       *pfname = filename;
767       if (strchr(filename, '*') && !fname_resolve(NULL, filename, pfname))
768       {
769         error("nema_hdrinfo_extract: no such file or filename", warning);
770         return(false);
771       }
772       
773       if ((store == NULL) || (store->to == NULL))
774         continue;
775       if ((fp = fopen_2(*pfname, "rb")) == NULL)
776         continue;
777       status = ((Bool (*) ())hdrextract_func)(fp, (Imrect *)(store->to));
778       fclose_2(fp, filename);
779       if (!status)
780         success = status;
781       store = store->next;
782     }
783 
784   return(success);
785 }
786 
787 
788 

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