imgfile.c
Go to the documentation of this file.
1 /******************************************************************************
2 
3  Copyright (c) 2003-2007 Turku PET Centre
4 
5  Library: imgfile.c
6  Description: I/O routines for IMG data.
7  Currently supported file formats:
8  ECAT 6.3 images and sinograms
9  ECAT 7.x 2D and 3D images (volumes) and sinograms
10  Analyze 7.5 images
11 
12  This library is free software; you can redistribute it and/or
13  modify it under the terms of the GNU Lesser General Public
14  License as published by the Free Software Foundation; either
15  version 2.1 of the License, or (at your option) any later version.
16 
17  This library is distributed in the hope that it will be useful,
18  but WITHOUT ANY WARRANTY; without even the implied warranty of
19  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
20  See the GNU Lesser General Public License for more details:
21  http://www.gnu.org/copyleft/lesser.html
22 
23  You should have received a copy of the GNU Lesser General Public License
24  along with this library/program; if not, write to the Free Software
25  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26 
27  Turku PET Centre, Turku, Finland, http://www.turkupetcentre.fi/
28 
29  Modification history:
30  2003-07-27 Vesa Oikonen
31  First created.
32  2003-08-11 VO
33  Corrected a bug in test print in imgRead().
34  2003-09-08 VO
35  Can now read&write 3D sinograms.
36  sampleDistance is now read/written from/to bin_size in ECAT7 mainheader.
37  imgWriteEcat7() renamed to imgWrite2DEcat7().
38  "New" function imgWriteEcat7() makes 3D format instead of 2D as before.
39  Applies the new IMG field _fileFormat in read/write.
40  2003-10-05 VO
41  Added support for Analyze 7.5 image format.
42  2003-10-07 VO
43  Setting Analyze db_name is corrected.
44  2003-10-08 VO
45  Only 7 first letters of image magic number is checked.
46  2003-10-21 VO
47  imgRead() and imgReadAnalyze() may work even if filename is given
48  with extension.
49  2003-11-04 VO
50  If ECAT 6.3 image mainheader and subheader contain the same calibration
51  factor, it is used only once.
52  Works with unit nCi/mL.
53  2003-11-10 VO
54  ECAT 7 3D scan subheader field x_resolution is used as sample distance.
55  2003-11-12 VO
56  Reading ECAT 6.3 image: pixels x size is read from subheader slice_width
57  if it is not found in main header plane_separation field.
58  2003-11-30 VO
59  For now, calls temp_roundf() instead of roundf().
60  2003-12-05 VO
61  Analyze files may be read&write without flipping in x,y,z-directions.
62  Function anaFlipping() is used to determine whether to flip or not.
63  2003-12-10 VO
64  Changes in 2003-11-10 and 2003-11-12 had been accidentally deleted
65  and are now returned.
66  2003-12-14 VO
67  Patient orientation is read&written in ECAT7 format.
68  Scanner (system_type) is read&written in ECAT formats.
69  2004-02-05 VO
70  anaFlipping() determines whether image is flipped in z-direction;
71  image is always flipped in x,y-directions.
72  2004-02-08 VO
73  Changes in imgSetEcat7MHeader():
74  -sets different magic number for sinograms and image volumes
75  -sets sw_version=72.
76  When writing ECAT 6.3 format: sets sw_version=2.
77  2004-02-22 VO
78  Analyze format: zoom factor is not written into/read from funused2,
79  because that field is used by SPM2 to another purpose.
80  2004-05-23 VO
81  ECAT7_3DSCANFIT is now supported as ECAT7_3DSCAN.
82  imgReadEcat7() reports that axial compression is not supported.
83  ECAT unit 1 is assumed to represent MBq/mL (IMG unit 13).
84  2004-05-24 VO
85  Pixel sizes are correctly converted from mm to cm when writing ECAT7.
86  Changes in unit conversions between ECAT and IMG.
87  2004-06-27 VO
88  Additional error message.
89  ecat63ReadAllToImg() and ecat63ReadPlaneToImg() do not even try to read
90  extra frames.
91  ECAT63_TEST changed to IMG_TEST.
92  Patient name is copied with strncpy in ecat63ReadPlaneToImg() and
93  ecat63ReadAllToImg().
94  2004-07-10 VO
95  Not so many test prints in reading ECAT 6.3 files.
96  2004-08-23 VO
97  MAX_STUDYNR_LEN applied where appropriate.
98  2004-09-20 VO
99  Doxygen style comments are corrected.
100  2004-09-24 VO
101  Better handling of ECAT7 calibration units.
102  E.g. imgUnitToEcat7() divided into imgUnitToEcat6() and imgUnitToEcat7().
103  2004-10-10 VO
104  imgSetEcat7MHeader() sets ECAT7 mainheader plane number to dimz, not 1
105  as before.
106  2004-10-13 VO
107  tm_isdst=-1 (unknown Daylight saving time) when appropriate.
108  2004-11-07
109  imgReadEcat7() reads files with axial compression, each slice as one plane.
110  2004-12-22 VO
111  Correction in imgGetEcat7MHeader(): calls imgUnitFromEcat7() instead of
112  imgUnitFromEcat().
113  2004-12-27 VO
114  imgUnitFromEcat(): IMG unit is set to unknown, not to MBq/mL, when ECAT
115  unit is unknown.
116  2005-03-03 Jarkko Johansson
117  initSIF changed to sifInit and readSIF to sifRead
118  2005-05-12 Calle Laakkonen
119  made image loading/saving functions filename argument const
120  2005-06-06 Calle Laakkonen
121  imgmsg is now shared by all functions.
122  2005-06-30 Harri Merisaari
123  fixed imgWrite(char* fname, IMG* img) to write also image structures
124  with '_fileFormat' field set as 'IMG_ANA_L'
125  2005-10-10 Calle Laakkonen
126  imgWriteAnalyze() now conserves memory by writing only 1 frame at a time.
127  2005-12-12 VO
128  imgReadAnalyze() sets img.isotopeHalflife, if isotope is found in SIF.
129  Corrected a setting of error message in imgRead().
130  2006-10-31 VO
131  Calibration unit functions moved to imgunit.c.
132  Return value of mktime() is checked.
133  Added "if(timezone == -7200) img.scanStart += 10800;"
134  into imgGetEcat7MHeader() as suggested by HM.
135  2007-01-31 VO
136  Corrected a bug in imgReadEcat7() with ECAT7 3D sinograms: all planes, not
137  only the first, are multiplied with dead-time correction factor.
138  Analyze 7.5 I/O functions separated into imgana.c.
139  ECAT 6.3 I/O functions separated into img_e63.c.
140  IMG status messages separated into imgmsg.h.
141  If valid study number is found in ECAT7 user_process_code,
142  then take it from there.
143  Patient_id and study_description are read and written in ECAT7.
144  Image x,y,z-resolutions are read/written in ECAT7.
145  Prompts and randoms (delayed) are read and written in ECAT7.
146  2007-02-27 VO
147  ECAT 7 I/O functions separated into img_e7.c.
148  2007-03-25 VO
149  Added functions imgFormatFromFName(), imgReadHeader(), imgReadNextFrame(),
150  imgReadFrame(), and imgWriteFrame().
151  2007-03-27 VO
152  imgFormatFromFName() identifies extension .polmap as polar map.
153  imgWrite() calls imgFormatFromFName() when necessary.
154  Read and write functions accept polar maps.
155  2007-17-07 HM
156  Changed void imgFormatFromFName(IMG *img, const char *fname) for ANSI
157 
158 ******************************************************************************/
159 #include <stdio.h>
160 #include <stdlib.h>
161 #include <unistd.h>
162 #include <math.h>
163 #include <string.h>
164 #include <time.h>
165 /*****************************************************************************/
166 #include "petc99.h"
167 #include "swap.h"
168 #include "halflife.h"
169 #include "substitutions.h"
170 /*****************************************************************************/
171 #include "include/img.h"
172 #include "include/ecat63.h"
173 #include "include/ecat7.h"
174 #include "include/analyze.h"
175 #include "include/imgmax.h"
176 #include "include/imgdecay.h"
177 #include "include/sif.h"
178 #include "include/imgfile.h"
179 /*****************************************************************************/
180 
181 /*****************************************************************************/
190 int imgRead(const char *fname, IMG *img) {
191  FILE *fp;
192  int ret;
193  ECAT7_mainheader ecat7_main_header;
194  ECAT63_mainheader ecat63_main_header;
195  char temp[FILENAME_MAX], *cptr;
196 
197  if(IMG_TEST) printf("imgRead(%s, *img)\n", fname);
198  /* Check the arguments */
199  if(fname==NULL) {img->statmsg=imgStatus(STATUS_FAULT); return(1);}
200  if(img==NULL || img->status!=IMG_STATUS_INITIALIZED) {
201  img->statmsg=imgStatus(STATUS_FAULT); return(2);}
202 
203  /* Check if we have Analyze file */
204  /* Check if filename was given accidentally with extension */
205  strcpy(temp, fname); cptr=strrchr(temp, '.');
206  if(cptr!=NULL && (strcmp(cptr, ".img")==0 || strcmp(cptr, ".hdr")==0))
207  *cptr=(char)0;
208  /* Are there the Analyze .img and .hdr file? */
209  if(anaExists(temp)!=0) {
210  /* Read Analyze image */
211  ret=imgReadAnalyze(temp, img);
212  if(ret==3 || ret==4 || ret==5 || ret==6) {
213  img->statmsg=imgStatus(STATUS_UNKNOWNFORMAT); return(4);
214  } else if(ret>0) {
215  img->statmsg=imgStatus(STATUS_NOFILE); return(4);
216  }
217  if(IMG_TEST) printf("%s identified as supported Analyze 7.5 format.\n",
218  fname);
220  return(0);
221  }
222 
223 
224  /* Check if we have an ECAT file */
225  /* Open file for read */
226  if((fp=fopen(fname, "rb")) == NULL) {
227  img->statmsg=imgStatus(STATUS_NOFILE); return(4);
228  }
229  /* Try to read ECAT 7.x main header */
230  ret=ecat7ReadMainheader(fp, &ecat7_main_header);
231  if(ret) {fclose(fp); img->statmsg=imgStatus(STATUS_UNKNOWNFORMAT); return(4);}
232  /* If header could be read, check for magic number */
233  if(strncmp(ecat7_main_header.magic_number, ECAT7V_MAGICNR, 7)==0) {
234  /* This is ECAT 7.x file */
235  /* Check if file type is supported */
236  if(imgEcat7Supported(&ecat7_main_header)==0) {
237  fclose(fp); img->statmsg=imgStatus(STATUS_UNSUPPORTED); return(5);
238  }
239  fclose(fp);
240  /* Read file */
241  if(IMG_TEST) printf("%s identified as supported ECAT 7.x %s format\n",
242  fname, ecat7filetype(ecat7_main_header.file_type));
243  ret=imgReadEcat7(fname, img);
244  if(ret) {if(IMG_TEST) printf("imgReadEcat7()=%d\n", ret); return(6);}
245  } else {
246  /* Check if file is in ECAT 6.3 format */
247  ret=ecat63ReadMainheader(fp, &ecat63_main_header);
248  fclose(fp);
249  if(ret==0) {
250  /* It seems to be ECAT 6.3, so read it */
251  if(IMG_TEST) printf("%s identified as supported ECAT 6.3 %s format\n",
252  fname, ecat7filetype(ecat63_main_header.file_type));
253  ret=ecat63ReadAllToImg(fname, img);
254  if(ret) {
255  if(IMG_TEST) fprintf(stderr, "ecat63ReaddAllToImg: %s\n", ecat63errmsg);
256  if(ret==6) img->statmsg=imgStatus(STATUS_MISSINGMATRIX);
258  return(6);
259  }
260  } else {img->statmsg=imgStatus(STATUS_UNKNOWNFORMAT); return(4);}
261  }
263  return(0);
264 }
265 /*****************************************************************************/
266 
267 /*****************************************************************************/
277 int imgWrite(const char *fname, IMG *img) {
278  int ret;
279 
280  if(IMG_TEST) printf("imgWrite(%s, *img)\n", fname);
281  /* Check the arguments */
282  if(fname==NULL) return(1);
283  if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) {
284  imgSetStatus(img, STATUS_FAULT); return(2);}
285  if(img->type!=IMG_TYPE_RAW &&
286  img->type!=IMG_TYPE_IMAGE &&
287  img->type!=IMG_TYPE_POLARMAP) {
288  imgSetStatus(img, STATUS_FAULT); return(2);}
289 
290  /* If _fileFormat is not defined, then determine it from the filename */
291  if(img->_fileFormat==IMG_UNKNOWN) imgFormatFromFName(img, fname);
292 
293  /* Write */
294  if(img->_fileFormat==IMG_E63) {
295  ret=ecat63WriteAllImg(fname, img);
296  switch(ret) {
297  case 0: break;
298  case 4: imgSetStatus(img, STATUS_NOMEMORY); break;
299  case 3: imgSetStatus(img, STATUS_NOWRITEPERM); break;
300  case 9: imgSetStatus(img, STATUS_DISKFULL); break;
301  default: imgSetStatus(img, STATUS_FAULT);
302  }
303  if(ret) return(7);
304  } else if(img->_fileFormat==IMG_ANA || img->_fileFormat==IMG_ANA_L) {
305  ret=imgWriteAnalyze(fname, img); if(ret) return(5);
306  } else if(img->_fileFormat==IMG_E7_2D) {
307  ret=imgWrite2DEcat7(fname, img); if(ret) return(5);
308  } else if(img->_fileFormat==IMG_POLARMAP) {
309  ret=imgWritePolarmap(fname, img); if(ret) return(5);
310  } else {
311  ret=imgWriteEcat7(fname, img); if(ret) return(5);
312  }
313  imgSetStatus(img, STATUS_OK);
314  return(0);
315 }
316 /*****************************************************************************/
317 
318 /*****************************************************************************/
329 int imgReadHeader(const char *fname, IMG *img) {
330  int ret;
331 
332 
333  if(IMG_TEST) printf("\nimgReadHeader(%s, *img)\n", fname);
334 
335  /* Check the arguments */
336  if(fname==NULL) return STATUS_FAULT;
337  if(img==NULL) return STATUS_FAULT;
338  if(img->status!=IMG_STATUS_INITIALIZED) return STATUS_FAULT;
339 
340  /* Check if we have Analyze database files */
341  if(anaDatabaseExists(fname, NULL, NULL, NULL)>0) { /* yes we have Analyze 7.5 */
342  /* Read Analyze header information */
343  ret=imgReadAnalyzeHeader(fname, img);
344  imgSetStatus(img, ret);
345  return(ret);
346  }
347 
348  /* Is this an ECAT7 file */
349  ret=imgReadEcat7Header(fname, img);
350  /* If main header was read but format was not identified as Ecat7,
351  it might be in Ecat6 format */
352  if(ret==STATUS_UNKNOWNFORMAT) {
353  /* Is this an ECAT6 file; check this as the last option, because
354  ECAT6 files don't contain any magic number */
355  ret=imgReadEcat63Header(fname, img);
356  if(ret) return STATUS_UNKNOWNFORMAT; /* don't know what this is */
357  /* Fine, it is Ecat 6.3 (or close enough) */
358  imgSetStatus(img, STATUS_OK);
359  return STATUS_OK;
360  }
361 
362  imgSetStatus(img, ret);
363  return ret;
364 }
365 /*****************************************************************************/
366 
367 /*****************************************************************************/
389 int imgReadFrame(const char *fname, int frame_to_read, IMG *img, int frame_index) {
390  IMG test_img;
391  int ret=0;
392 
393  if(IMG_TEST) printf("\nimgReadFrame(%s, %d, *img, %d)\n",
394  fname, frame_to_read, frame_index);
395 
396  /*
397  * Check the input
398  */
399  if(fname==NULL) return STATUS_FAULT;
400  if(img==NULL) return STATUS_FAULT;
402  return STATUS_FAULT;
403  if(frame_to_read<1) return STATUS_FAULT;
404  if(frame_index<0) return STATUS_FAULT;
405  /* if frame_index>0, then there must be sufficient memory allocated for it */
406  if(frame_index>0) {
407  if(img->status!=IMG_STATUS_OCCUPIED) return STATUS_FAULT;
408  if(frame_index>img->dimt-1) return STATUS_FAULT;
409  }
410 
411  /*
412  * If IMG is preallocated, check that fundamental header information
413  * is compatible with old and new contents.
414  * If not allocated, then read the header contents, and allocate it
415  */
416  imgInit(&test_img);
417  if(img->status==IMG_STATUS_OCCUPIED) {
418  ret=imgReadHeader(fname, &test_img); imgSetStatus(&test_img, ret);
419  if(IMG_TEST>1) printf("imgReadHeader() return message := %s\n", test_img.statmsg);
420  if(ret) return(ret);
421  if(IMG_TEST>3) imgInfo(&test_img);
422  /* Test that file format and type are the same */
423  ret=0;
424  if(img->type!=test_img.type) ret++;
425  if(img->_fileFormat!=test_img._fileFormat) ret++;
426  /* Test that x, y, and z dimensions are the same */
427  if(img->dimx!=test_img.dimx) ret++;
428  if(img->dimy!=test_img.dimy) ret++;
429  if(img->dimz!=test_img.dimz) ret++;
430  imgEmpty(&test_img); if(ret>0) return STATUS_INVALIDHEADER;
431  } else {
432  ret=imgReadHeader(fname, img); imgSetStatus(img, ret);
433  if(IMG_TEST>1) printf("imgReadHeader() return message := %s\n", img->statmsg);
434  if(ret) return(ret);
435  if(IMG_TEST>3) imgInfo(img);
436  /* Allocate memory for one frame */
437  img->dimt=1;
438  ret=imgAllocate(img, img->dimz, img->dimy, img->dimx, img->dimt);
439  if(ret) return STATUS_NOMEMORY;
440  }
441 
442  /*
443  * Read the frame data and corresponding information like frame time
444  * if available
445  */
446  switch(img->_fileFormat) {
447  case IMG_E7:
448  case IMG_E7_2D:
449  case IMG_POLARMAP:
450  ret=imgReadEcat7Frame(fname, frame_to_read, img, frame_index);
451  if(IMG_TEST>1) printf("imgReadEcat7Frame() return value := %d\n", ret);
452  break;
453  case IMG_E63:
454  ret=imgReadEcat63Frame(fname, frame_to_read, img, frame_index);
455  if(IMG_TEST>1) printf("imgReadEcat63Frame() return value := %d\n", ret);
456  break;
457  case IMG_ANA:
458  case IMG_ANA_L:
459  ret=imgReadAnalyzeFrame(fname, frame_to_read, img, frame_index);
460  if(IMG_TEST>1) printf("imgReadAnalyzeFrame() return value := %d\n", ret);
461  break;
462  default:
463  ret=STATUS_UNSUPPORTED;
464  }
465  imgSetStatus(img, ret);
466  return ret;
467 }
468 /*****************************************************************************/
469 
470 /*****************************************************************************/
493 int imgWriteFrame(const char *fname, int frame_to_write, IMG *img, int frame_index) {
494  int ret=0;
495 
496  if(IMG_TEST) printf("\nimgWriteFrame(%s, %d, *img, %d)\n",
497  fname, frame_to_write, frame_index);
498 
499  /*
500  * Check the input
501  */
502  if(fname==NULL) return STATUS_FAULT;
503  if(img==NULL) return STATUS_FAULT;
504  if(img->status!=IMG_STATUS_OCCUPIED) return STATUS_FAULT;
505  if(frame_to_write<0) return STATUS_FAULT;
506  if(frame_index<0 || frame_index>=img->dimt) return STATUS_FAULT;
507 
508 
509  /*
510  * Call separate function for each supported file format
511  */
512  imgFormatFromFName(img, fname);
513  switch(img->_fileFormat) {
514  case IMG_E7:
515  case IMG_E7_2D:
516  case IMG_POLARMAP:
517  ret=imgWriteEcat7Frame(fname, frame_to_write, img, frame_index);
518  break;
519  case IMG_E63:
520  ret=imgWriteEcat63Frame(fname, frame_to_write, img, frame_index);
521  break;
522  case IMG_ANA:
523  case IMG_ANA_L:
524  ret=STATUS_UNSUPPORTED;
525  /* Not supported because would require global min&max values
526  * if saved in short ints which is now the only possibility
527  * ret=imgWriteAnaFrame(fname, frame_to_write, img, frame_index);
528  */
529  break;
530  default:
531  ret=STATUS_UNSUPPORTED;
532  }
533  imgSetStatus(img, ret);
534  return ret;
535 }
536 /*****************************************************************************/
537 
538 /*****************************************************************************/
547 void imgFormatFromFName(IMG *img, const char *fname) {
548  char *cptr=NULL;
549 
550  if(img->_fileFormat!=IMG_UNKNOWN && img->_fileFormat>0) return;
551  img->_fileFormat=IMG_E7; /* default */
552  /* get extension */
553  cptr=strrchr(fname, '.'); if(cptr!=NULL) cptr++;
554  if(cptr!=NULL) {
555  if(strcasecmp(cptr, "hdr")==0) { img->_fileFormat=IMG_ANA; return;}
556  if(strcasecmp(cptr, "polmap")==0) { img->_fileFormat=IMG_POLARMAP; return;}
557  if(strcasecmp(cptr, "img")==0 ||
558  strcasecmp(cptr, "scn")==0 ||
559  strcasecmp(cptr, "nrm")==0 ||
560  strcasecmp(cptr, "atn")==0) {
561  img->_fileFormat=IMG_E63; return;}
562  } else { /* no extension at all */
563  img->_fileFormat=IMG_ANA;
564  }
565 }
566 /*****************************************************************************/
567 
568 /*****************************************************************************/