imgfile.c
Go to the documentation of this file.
00001 /******************************************************************************
00002 
00003   Copyright (c) 2003-2007 Turku PET Centre
00004 
00005   Library:     imgfile.c
00006   Description: I/O routines for IMG data.
00007                Currently supported file formats:
00008                ECAT 6.3 images and sinograms
00009                ECAT 7.x 2D and 3D images (volumes) and sinograms
00010                Analyze 7.5 images
00011 
00012   This library is free software; you can redistribute it and/or
00013   modify it under the terms of the GNU Lesser General Public
00014   License as published by the Free Software Foundation; either
00015   version 2.1 of the License, or (at your option) any later version.
00016 
00017   This library is distributed in the hope that it will be useful,
00018   but WITHOUT ANY WARRANTY; without even the implied warranty of
00019   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
00020   See the GNU Lesser General Public License for more details:
00021   http://www.gnu.org/copyleft/lesser.html
00022 
00023   You should have received a copy of the GNU Lesser General Public License
00024   along with this library/program; if not, write to the Free Software
00025   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 
00026 
00027   Turku PET Centre, Turku, Finland, http://www.turkupetcentre.fi/
00028 
00029   Modification history:
00030   2003-07-27 Vesa Oikonen
00031     First created.
00032   2003-08-11 VO
00033     Corrected a bug in test print in imgRead().
00034   2003-09-08 VO
00035     Can now read&write 3D sinograms.
00036     sampleDistance is now read/written from/to bin_size in ECAT7 mainheader.
00037     imgWriteEcat7() renamed to imgWrite2DEcat7().
00038     "New" function imgWriteEcat7() makes 3D format instead of 2D as before.
00039     Applies the new IMG field _fileFormat in read/write.
00040   2003-10-05 VO
00041     Added support for Analyze 7.5 image format.
00042   2003-10-07 VO
00043     Setting Analyze db_name is corrected.
00044   2003-10-08 VO
00045     Only 7 first letters of image magic number is checked.
00046   2003-10-21 VO
00047     imgRead() and imgReadAnalyze() may work even if filename is given
00048     with extension.
00049   2003-11-04 VO
00050     If ECAT 6.3 image mainheader and subheader contain the same calibration
00051     factor, it is used only once.
00052     Works with unit nCi/mL.
00053   2003-11-10 VO
00054     ECAT 7 3D scan subheader field x_resolution is used as sample distance.
00055   2003-11-12 VO
00056     Reading ECAT 6.3 image: pixels x size is read from subheader slice_width
00057     if it is not found in main header plane_separation field.
00058   2003-11-30 VO
00059     For now, calls temp_roundf() instead of roundf().
00060   2003-12-05 VO
00061     Analyze files may be read&write without flipping in x,y,z-directions.
00062     Function anaFlipping() is used to determine whether to flip or not.
00063   2003-12-10 VO
00064     Changes in 2003-11-10 and 2003-11-12 had been accidentally deleted
00065     and are now returned.
00066   2003-12-14 VO
00067     Patient orientation is read&written in ECAT7 format.
00068     Scanner (system_type) is read&written in ECAT formats.
00069   2004-02-05 VO
00070     anaFlipping() determines whether image is flipped in z-direction;
00071     image is always flipped in x,y-directions.
00072   2004-02-08 VO
00073     Changes in imgSetEcat7MHeader():
00074     -sets different magic number for sinograms and image volumes
00075     -sets sw_version=72.
00076     When writing ECAT 6.3 format: sets sw_version=2.
00077   2004-02-22 VO
00078     Analyze format: zoom factor is not written into/read from funused2,
00079     because that field is used by SPM2 to another purpose.
00080   2004-05-23 VO
00081     ECAT7_3DSCANFIT is now supported as ECAT7_3DSCAN.
00082     imgReadEcat7() reports that axial compression is not supported.
00083     ECAT unit 1 is assumed to represent MBq/mL (IMG unit 13).
00084   2004-05-24 VO
00085     Pixel sizes are correctly converted from mm to cm when writing ECAT7.
00086     Changes in unit conversions between ECAT and IMG.
00087   2004-06-27 VO
00088     Additional error message.
00089     ecat63ReadAllToImg() and ecat63ReadPlaneToImg() do not even try to read
00090     extra frames.
00091     ECAT63_TEST changed to IMG_TEST.
00092     Patient name is copied with strncpy in ecat63ReadPlaneToImg() and
00093     ecat63ReadAllToImg().
00094   2004-07-10 VO
00095     Not so many test prints in reading ECAT 6.3 files.
00096   2004-08-23 VO
00097     MAX_STUDYNR_LEN applied where appropriate.
00098   2004-09-20 VO
00099     Doxygen style comments are corrected.
00100   2004-09-24 VO
00101     Better handling of ECAT7 calibration units.
00102     E.g. imgUnitToEcat7() divided into imgUnitToEcat6() and imgUnitToEcat7().
00103   2004-10-10 VO
00104     imgSetEcat7MHeader() sets ECAT7 mainheader plane number to dimz, not 1
00105     as before.
00106   2004-10-13 VO
00107     tm_isdst=-1 (unknown Daylight saving time) when appropriate.
00108   2004-11-07
00109     imgReadEcat7() reads files with axial compression, each slice as one plane.
00110   2004-12-22 VO
00111     Correction in imgGetEcat7MHeader(): calls imgUnitFromEcat7() instead of
00112     imgUnitFromEcat().
00113   2004-12-27 VO
00114     imgUnitFromEcat(): IMG unit is set to unknown, not to MBq/mL, when ECAT
00115     unit is unknown.
00116   2005-03-03 Jarkko Johansson
00117     initSIF changed to sifInit and readSIF to sifRead
00118   2005-05-12 Calle Laakkonen
00119     made image loading/saving functions filename argument const
00120   2005-06-06 Calle Laakkonen
00121     imgmsg is now shared by all functions.
00122   2005-06-30 Harri Merisaari
00123     fixed imgWrite(char* fname, IMG* img) to write also image structures 
00124     with '_fileFormat' field set as 'IMG_ANA_L' 
00125   2005-10-10 Calle Laakkonen
00126     imgWriteAnalyze() now conserves memory by writing only 1 frame at a time.
00127   2005-12-12 VO
00128     imgReadAnalyze() sets img.isotopeHalflife, if isotope is found in SIF.
00129     Corrected a setting of error message in imgRead().
00130   2006-10-31 VO
00131     Calibration unit functions moved to imgunit.c.
00132     Return value of mktime() is checked.
00133     Added "if(timezone == -7200) img.scanStart += 10800;"
00134     into imgGetEcat7MHeader() as suggested by HM.
00135   2007-01-31 VO
00136     Corrected a bug in imgReadEcat7() with ECAT7 3D sinograms: all planes, not
00137     only the first, are multiplied with dead-time correction factor.
00138     Analyze 7.5 I/O functions separated into imgana.c.
00139     ECAT 6.3 I/O functions separated into img_e63.c.
00140     IMG status messages separated into imgmsg.h.
00141     If valid study number is found in ECAT7 user_process_code, 
00142     then take it from there.  
00143     Patient_id and study_description are read and written in ECAT7.
00144     Image x,y,z-resolutions are read/written in ECAT7. 
00145     Prompts and randoms (delayed) are read and written in ECAT7.
00146   2007-02-27 VO
00147     ECAT 7 I/O functions separated into img_e7.c.
00148   2007-03-25 VO
00149     Added functions imgFormatFromFName(), imgReadHeader(), imgReadNextFrame(),
00150     imgReadFrame(), and imgWriteFrame().
00151   2007-03-27 VO
00152     imgFormatFromFName() identifies extension .polmap as polar map.
00153     imgWrite() calls imgFormatFromFName() when necessary.
00154     Read and write functions accept polar maps.
00155   2007-17-07 HM
00156     Changed void imgFormatFromFName(IMG *img, const char *fname) for ANSI
00157 
00158 ******************************************************************************/
00159 #include <stdio.h>
00160 #include <stdlib.h>
00161 #include <unistd.h>
00162 #include <math.h>
00163 #include <string.h>
00164 #include <time.h>
00165 /*****************************************************************************/
00166 #include "petc99.h"
00167 #include "swap.h"
00168 #include "halflife.h"
00169 #include "substitutions.h"
00170 /*****************************************************************************/
00171 #include "include/img.h"
00172 #include "include/ecat63.h"
00173 #include "include/ecat7.h"
00174 #include "include/analyze.h"
00175 #include "include/imgmax.h"
00176 #include "include/imgdecay.h"
00177 #include "include/sif.h"
00178 #include "include/imgfile.h"
00179 /*****************************************************************************/
00180 
00181 /*****************************************************************************/
00190 int imgRead(const char *fname, IMG *img) {
00191   FILE *fp;
00192   int ret;
00193   ECAT7_mainheader ecat7_main_header;
00194   ECAT63_mainheader ecat63_main_header;
00195   char temp[FILENAME_MAX], *cptr;
00196 
00197   if(IMG_TEST) printf("imgRead(%s, *img)\n", fname);
00198   /* Check the arguments */
00199   if(fname==NULL) {img->statmsg=imgStatus(STATUS_FAULT); return(1);}
00200   if(img==NULL || img->status!=IMG_STATUS_INITIALIZED) {
00201     img->statmsg=imgStatus(STATUS_FAULT); return(2);}
00202 
00203   /* Check if we have Analyze file */
00204   /* Check if filename was given accidentally with extension */
00205   strcpy(temp, fname); cptr=strrchr(temp, '.');
00206   if(cptr!=NULL && (strcmp(cptr, ".img")==0 || strcmp(cptr, ".hdr")==0))
00207     *cptr=(char)0;
00208   /* Are there the Analyze .img and .hdr file? */
00209   if(anaExists(temp)!=0) {
00210     /* Read Analyze image */
00211     ret=imgReadAnalyze(temp, img);
00212     if(ret==3 || ret==4 || ret==5 || ret==6) {
00213       img->statmsg=imgStatus(STATUS_UNKNOWNFORMAT); return(4);
00214     } else if(ret>0) {
00215       img->statmsg=imgStatus(STATUS_NOFILE); return(4);
00216     }
00217     if(IMG_TEST) printf("%s identified as supported Analyze 7.5 format.\n",
00218       fname);
00219     img->statmsg=imgStatus(STATUS_OK);
00220     return(0);
00221   }
00222 
00223 
00224   /* Check if we have an ECAT file */
00225   /* Open file for read */
00226   if((fp=fopen(fname, "rb")) == NULL) {
00227     img->statmsg=imgStatus(STATUS_NOFILE); return(4);
00228   }
00229   /* Try to read ECAT 7.x main header */
00230   ret=ecat7ReadMainheader(fp, &ecat7_main_header);
00231   if(ret) {fclose(fp); img->statmsg=imgStatus(STATUS_UNKNOWNFORMAT); return(4);}
00232   /* If header could be read, check for magic number */
00233   if(strncmp(ecat7_main_header.magic_number, ECAT7V_MAGICNR, 7)==0) {
00234     /* This is ECAT 7.x file */
00235     /* Check if file type is supported */
00236     if(imgEcat7Supported(&ecat7_main_header)==0) {
00237       fclose(fp); img->statmsg=imgStatus(STATUS_UNSUPPORTED); return(5);
00238     }
00239     fclose(fp);
00240     /* Read file */
00241     if(IMG_TEST) printf("%s identified as supported ECAT 7.x %s format\n",
00242       fname, ecat7filetype(ecat7_main_header.file_type));
00243     ret=imgReadEcat7(fname, img);
00244     if(ret) {if(IMG_TEST) printf("imgReadEcat7()=%d\n", ret); return(6);}
00245   } else {
00246     /* Check if file is in ECAT 6.3 format */
00247     ret=ecat63ReadMainheader(fp, &ecat63_main_header);
00248     fclose(fp);
00249     if(ret==0) {
00250       /* It seems to be ECAT 6.3, so read it */
00251       if(IMG_TEST) printf("%s identified as supported ECAT 6.3 %s format\n",
00252         fname, ecat7filetype(ecat63_main_header.file_type));
00253       ret=ecat63ReadAllToImg(fname, img);
00254       if(ret) {
00255         if(IMG_TEST) fprintf(stderr, "ecat63ReaddAllToImg: %s\n", ecat63errmsg);
00256         if(ret==6) img->statmsg=imgStatus(STATUS_MISSINGMATRIX);
00257   else img->statmsg=imgStatus(STATUS_UNSUPPORTED);
00258         return(6);
00259       }
00260     } else {img->statmsg=imgStatus(STATUS_UNKNOWNFORMAT); return(4);}
00261   }
00262   img->statmsg=imgStatus(STATUS_OK);
00263   return(0);
00264 }
00265 /*****************************************************************************/
00266 
00267 /*****************************************************************************/
00277 int imgWrite(const char *fname, IMG *img) {
00278   int ret;
00279 
00280   if(IMG_TEST) printf("imgWrite(%s, *img)\n", fname);
00281   /* Check the arguments */
00282   if(fname==NULL) return(1);
00283   if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) {
00284     imgSetStatus(img, STATUS_FAULT); return(2);}
00285   if(img->type!=IMG_TYPE_RAW &&
00286      img->type!=IMG_TYPE_IMAGE &&
00287      img->type!=IMG_TYPE_POLARMAP) {
00288     imgSetStatus(img, STATUS_FAULT); return(2);}
00289 
00290   /* If _fileFormat is not defined, then determine it from the filename */
00291   if(img->_fileFormat==IMG_UNKNOWN) imgFormatFromFName(img, fname);
00292 
00293   /* Write */
00294   if(img->_fileFormat==IMG_E63) {
00295     ret=ecat63WriteAllImg(fname, img);
00296     switch(ret) {
00297       case 0: break;
00298       case 4: imgSetStatus(img, STATUS_NOMEMORY); break;
00299       case 3: imgSetStatus(img, STATUS_NOWRITEPERM); break;
00300       case 9: imgSetStatus(img, STATUS_DISKFULL); break;
00301       default: imgSetStatus(img, STATUS_FAULT);
00302     }
00303     if(ret) return(7);
00304   } else if(img->_fileFormat==IMG_ANA || img->_fileFormat==IMG_ANA_L) {
00305     ret=imgWriteAnalyze(fname, img); if(ret) return(5);
00306   } else if(img->_fileFormat==IMG_E7_2D) {
00307     ret=imgWrite2DEcat7(fname, img); if(ret) return(5);
00308   } else if(img->_fileFormat==IMG_POLARMAP) {
00309     ret=imgWritePolarmap(fname, img); if(ret) return(5);
00310   } else {
00311     ret=imgWriteEcat7(fname, img); if(ret) return(5);
00312   }
00313   imgSetStatus(img, STATUS_OK);
00314   return(0);
00315 }
00316 /*****************************************************************************/
00317 
00318 /*****************************************************************************/
00329 int imgReadHeader(const char *fname, IMG *img) {
00330   int ret;
00331 
00332 
00333   if(IMG_TEST) printf("\nimgReadHeader(%s, *img)\n", fname);
00334 
00335   /* Check the arguments */
00336   if(fname==NULL) return STATUS_FAULT;
00337   if(img==NULL) return STATUS_FAULT;
00338   if(img->status!=IMG_STATUS_INITIALIZED) return STATUS_FAULT;
00339 
00340   /* Check if we have Analyze database files */
00341   if(anaDatabaseExists(fname, NULL, NULL, NULL)>0) { /* yes we have Analyze 7.5 */
00342     /* Read Analyze header information */
00343     ret=imgReadAnalyzeHeader(fname, img);
00344     imgSetStatus(img, ret);
00345     return(ret);
00346   }
00347 
00348   /* Is this an ECAT7 file */
00349   ret=imgReadEcat7Header(fname, img);
00350   /* If main header was read but format was not identified as Ecat7, 
00351      it might be in Ecat6 format */
00352   if(ret==STATUS_UNKNOWNFORMAT) {
00353     /* Is this an ECAT6 file; check this as the last option, because
00354        ECAT6 files don't contain any magic number */
00355     ret=imgReadEcat63Header(fname, img);
00356     if(ret) return STATUS_UNKNOWNFORMAT; /* don't know what this is */
00357     /* Fine, it is Ecat 6.3 (or close enough) */
00358     imgSetStatus(img, STATUS_OK);
00359     return STATUS_OK;
00360   }
00361   
00362   imgSetStatus(img, ret);
00363   return ret;
00364 }
00365 /*****************************************************************************/
00366 
00367 /*****************************************************************************/
00389 int imgReadFrame(const char *fname, int frame_to_read, IMG *img, int frame_index) {
00390   IMG test_img;
00391   int ret=0;
00392 
00393   if(IMG_TEST) printf("\nimgReadFrame(%s, %d, *img, %d)\n",
00394     fname, frame_to_read, frame_index);
00395     
00396   /*
00397    *  Check the input 
00398    */
00399   if(fname==NULL) return STATUS_FAULT;
00400   if(img==NULL) return STATUS_FAULT;
00401   if(img->status!=IMG_STATUS_INITIALIZED && img->status!=IMG_STATUS_OCCUPIED)
00402     return STATUS_FAULT;
00403   if(frame_to_read<1) return STATUS_FAULT;
00404   if(frame_index<0) return STATUS_FAULT;
00405   /* if frame_index>0, then there must be sufficient memory allocated for it */
00406   if(frame_index>0) {
00407     if(img->status!=IMG_STATUS_OCCUPIED) return STATUS_FAULT;
00408     if(frame_index>img->dimt-1) return STATUS_FAULT;
00409   }
00410 
00411   /*
00412    *  If IMG is preallocated, check that fundamental header information
00413    *  is compatible with old and new contents.
00414    *  If not allocated, then read the header contents, and allocate it
00415    */
00416   imgInit(&test_img);
00417   if(img->status==IMG_STATUS_OCCUPIED) {
00418     ret=imgReadHeader(fname, &test_img); imgSetStatus(&test_img, ret);
00419     if(IMG_TEST>1) printf("imgReadHeader() return message := %s\n", test_img.statmsg);
00420     if(ret) return(ret);
00421     if(IMG_TEST>3) imgInfo(&test_img);
00422     /* Test that file format and type are the same */
00423     ret=0;
00424     if(img->type!=test_img.type) ret++;
00425     if(img->_fileFormat!=test_img._fileFormat) ret++;
00426     /* Test that x, y, and z dimensions are the same */
00427     if(img->dimx!=test_img.dimx) ret++;
00428     if(img->dimy!=test_img.dimy) ret++;
00429     if(img->dimz!=test_img.dimz) ret++; 
00430     imgEmpty(&test_img); if(ret>0) return STATUS_INVALIDHEADER;
00431   } else {
00432     ret=imgReadHeader(fname, img); imgSetStatus(img, ret);
00433     if(IMG_TEST>1) printf("imgReadHeader() return message := %s\n", img->statmsg);
00434     if(ret) return(ret);
00435     if(IMG_TEST>3) imgInfo(img);
00436     /* Allocate memory for one frame */
00437     img->dimt=1;
00438     ret=imgAllocate(img, img->dimz, img->dimy, img->dimx, img->dimt);
00439     if(ret) return STATUS_NOMEMORY;
00440   }
00441 
00442   /*
00443    *  Read the frame data and corresponding information like frame time
00444    *  if available   
00445    */
00446   switch(img->_fileFormat) {
00447     case IMG_E7:
00448     case IMG_E7_2D:
00449     case IMG_POLARMAP:
00450       ret=imgReadEcat7Frame(fname, frame_to_read, img, frame_index);
00451       if(IMG_TEST>1) printf("imgReadEcat7Frame() return value := %d\n", ret);
00452       break;
00453     case IMG_E63:
00454       ret=imgReadEcat63Frame(fname, frame_to_read, img, frame_index);
00455       if(IMG_TEST>1) printf("imgReadEcat63Frame() return value := %d\n", ret);
00456       break;
00457     case IMG_ANA:
00458     case IMG_ANA_L:
00459       ret=imgReadAnalyzeFrame(fname, frame_to_read, img, frame_index);
00460       if(IMG_TEST>1) printf("imgReadAnalyzeFrame() return value := %d\n", ret);
00461       break;
00462     default:
00463       ret=STATUS_UNSUPPORTED;
00464   }
00465   imgSetStatus(img, ret);
00466   return ret;
00467 }
00468 /*****************************************************************************/
00469 
00470 /*****************************************************************************/
00493 int imgWriteFrame(const char *fname, int frame_to_write, IMG *img, int frame_index) {
00494   int ret=0;
00495 
00496   if(IMG_TEST) printf("\nimgWriteFrame(%s, %d, *img, %d)\n",
00497     fname, frame_to_write, frame_index);
00498     
00499   /*
00500    *  Check the input 
00501    */
00502   if(fname==NULL) return STATUS_FAULT;
00503   if(img==NULL) return STATUS_FAULT;
00504   if(img->status!=IMG_STATUS_OCCUPIED) return STATUS_FAULT;
00505   if(frame_to_write<0) return STATUS_FAULT;
00506   if(frame_index<0 || frame_index>=img->dimt) return STATUS_FAULT;
00507 
00508 
00509   /*
00510    *  Call separate function for each supported file format
00511    */
00512   imgFormatFromFName(img, fname);
00513   switch(img->_fileFormat) {
00514     case IMG_E7:
00515     case IMG_E7_2D:
00516     case IMG_POLARMAP:
00517       ret=imgWriteEcat7Frame(fname, frame_to_write, img, frame_index);
00518       break;
00519     case IMG_E63:
00520       ret=imgWriteEcat63Frame(fname, frame_to_write, img, frame_index);
00521       break;
00522     case IMG_ANA:
00523     case IMG_ANA_L:
00524       ret=STATUS_UNSUPPORTED;
00525       /* Not supported because would require global min&max values
00526        * if saved in short ints which is now the only possibility
00527        * ret=imgWriteAnaFrame(fname, frame_to_write, img, frame_index);
00528        */
00529       break;
00530     default:
00531       ret=STATUS_UNSUPPORTED;
00532   }
00533   imgSetStatus(img, ret);
00534   return ret;
00535 }
00536 /*****************************************************************************/
00537 
00538 /*****************************************************************************/
00547 void imgFormatFromFName(IMG *img, const char *fname) {
00548   char *cptr=NULL;
00549 
00550   if(img->_fileFormat!=IMG_UNKNOWN && img->_fileFormat>0) return;
00551   img->_fileFormat=IMG_E7; /* default */
00552   /* get extension */
00553   cptr=strrchr(fname, '.'); if(cptr!=NULL) cptr++;
00554   if(cptr!=NULL) {
00555     if(strcasecmp(cptr, "hdr")==0) { img->_fileFormat=IMG_ANA; return;}
00556     if(strcasecmp(cptr, "polmap")==0) { img->_fileFormat=IMG_POLARMAP; return;}
00557     if(strcasecmp(cptr, "img")==0 ||
00558        strcasecmp(cptr, "scn")==0 ||
00559        strcasecmp(cptr, "nrm")==0 ||
00560        strcasecmp(cptr, "atn")==0) {
00561       img->_fileFormat=IMG_E63; return;}
00562   } else { /* no extension at all */
00563     img->_fileFormat=IMG_ANA;
00564   }
00565 }
00566 /*****************************************************************************/
00567 
00568 /*****************************************************************************/