img_ana.c
Go to the documentation of this file.
00001 /******************************************************************************
00002 
00003   Copyright (c) 2007,2008,2009 Turku PET Centre
00004 
00005   Library:     img_ana.c
00006   Description: I/O routines for IMG data from/to Analyze 7.5 format.
00007 
00008   This library is free software; you can redistribute it and/or
00009   modify it under the terms of the GNU Lesser General Public
00010   License as published by the Free Software Foundation; either
00011   version 2.1 of the License, or (at your option) any later version.
00012 
00013   This library is distributed in the hope that it will be useful,
00014   but WITHOUT ANY WARRANTY; without even the implied warranty of
00015   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
00016   See the GNU Lesser General Public License for more details:
00017   http://www.gnu.org/copyleft/lesser.html
00018 
00019   You should have received a copy of the GNU Lesser General Public License
00020   along with this library/program; if not, write to the Free Software
00021   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 
00022 
00023   Turku PET Centre, Turku, Finland, http://www.turkupetcentre.fi/
00024 
00025   Modification history:
00026   2007-01-31 Vesa Oikonen
00027     Functions moved from imgfile.c.
00028     Prompts and randoms are read from SIF data to IMG when reading Analyze.
00029   2007-02-27 VO
00030     Bug correction in imgWriteAnalyze(): fp was not closed in all errors.
00031   2007-03-25 VO
00032     Added functions imgReadAnalyzeHeader(), imgGetAnalyzeHeader(),
00033     imgSetAnalyzeHeader(), imgReadAnalyzeFrame(), imgReadAnalyzeFirstFrame(),
00034     and imgWriteAnalyzeFrame().
00035   2007-17-07 Harri Merisaari
00036     Modified for optional ANSi compatibility    
00037   2007-09-10 VO
00038     Return value of localtime() is checked.
00039   2008-07-07 VO
00040     If information on decay correction is not found in Analyze header, then
00041     it is assumed that image data is corrected for decay; previously assumed
00042     that image was NOT corrected for decay.
00043   2009-12-10 VO
00044     strcpy() replaced with strncpy() in filling header fields.
00045 
00046 
00047 ******************************************************************************/
00048 #include <stdio.h>
00049 #include <stdlib.h>
00050 #include <unistd.h>
00051 #include <math.h>
00052 #include <string.h>
00053 #include <time.h>
00054 /*****************************************************************************/
00055 #include "petc99.h"
00056 #include "swap.h"
00057 #include "halflife.h"
00058 /*****************************************************************************/
00059 #include "include/img.h"
00060 #include "include/analyze.h"
00061 #include "include/imgmax.h"
00062 #include "include/imgdecay.h"
00063 #include "include/sif.h"
00064 #include "include/imgfile.h"
00065 /*****************************************************************************/
00066 
00067 /*****************************************************************************/
00083 int imgReadAnalyze(const char *dbname, IMG *img) {
00084   FILE *fp;
00085   int ret, fi, pi, xi, yi;
00086   float *fdata=NULL, *fptr;
00087   ANALYZE_DSR dsr;
00088   char datfile[FILENAME_MAX], hdrfile[FILENAME_MAX], siffile[FILENAME_MAX];
00089   char buf[128], *cptr;
00090   int dimNr, dimx, dimy, dimz=1, dimt=1, pxlNr=0;
00091   SIF sif;
00092   struct tm *st;
00093 
00094 
00095   if(IMG_TEST) printf("imgReadAnalyze(%s, *img)\n", dbname);
00096 
00097   /* Check the arguments */
00098   imgSetStatus(img, STATUS_OK);
00099   if(img==NULL || img->status!=IMG_STATUS_INITIALIZED) {
00100     imgSetStatus(img, STATUS_FAULT); return(2);}
00101   if(dbname==NULL || !dbname[0]) {imgSetStatus(img, STATUS_FAULT); return(1);}
00102   
00103   /* Make the image and header filenames */
00104   if(anaExists(dbname)==0) {
00105     /* Check if filename was given accidentally with extension */
00106     strcpy(datfile, dbname); cptr=strrchr(datfile, '.');
00107     if(cptr!=NULL && (strcmp(cptr, ".img")==0 || strcmp(cptr, ".hdr")==0)) {
00108       *cptr=(char)0; strcpy(hdrfile, datfile);
00109       if(anaExists(datfile)==0) { /* still not found */
00110         imgSetStatus(img, STATUS_NOHEADERFILE); return(3);}
00111       strcat(datfile, ".img"); strcat(hdrfile, ".hdr");
00112     } else {
00113       imgSetStatus(img, STATUS_NOHEADERFILE); return(3);
00114     }
00115   } else {
00116     /* Database name was given and img and hdr files were found */
00117     strcpy(datfile, dbname); strcat(datfile, ".img");
00118     strcpy(hdrfile, dbname); strcat(hdrfile, ".hdr");
00119   }
00120 
00121   /* Read Analyze header file */
00122   ret=anaReadHeader(hdrfile, &dsr);
00123   if(ret) {
00124     if(ret==1) imgSetStatus(img, STATUS_FAULT);
00125     else if(ret==2) imgSetStatus(img, STATUS_NOHEADERFILE);
00126     else imgSetStatus(img, STATUS_UNSUPPORTED);
00127     return(3);
00128   }
00129   if(IMG_TEST) anaPrintHeader(&dsr, stdout);
00130 
00131   /* Open image datafile */
00132   if(IMG_TEST) fprintf(stdout, "reading image data %s\n", datfile);
00133   if((fp=fopen(datfile, "rb")) == NULL) {imgSetStatus(img, STATUS_NOIMGDATA); return(5);}
00134 
00135   /* Prepare IMG for Analyze image */
00136   /* Get the image dimensions from header */
00137   dimNr=dsr.dime.dim[0];
00138   if(dimNr<2) {fclose(fp); imgSetStatus(img, STATUS_INVALIDHEADER); return(4);}
00139   dimx=dsr.dime.dim[1]; dimy=dsr.dime.dim[2];
00140   if(dimNr>2) {dimz=dsr.dime.dim[3]; if(dimNr>3) dimt=dsr.dime.dim[4];}
00141   pxlNr=dimx*dimy*dimz;
00142   if(pxlNr<1) {fclose(fp); imgSetStatus(img, STATUS_INVALIDHEADER); return(4);}
00143   /* Allocate memory for IMG */
00144   ret=imgAllocate(img, dimz, dimy, dimx, dimt);
00145   if(ret) {fclose(fp); imgSetStatus(img, STATUS_NOMEMORY); return(11);}
00146   /* Copy information from Analyze header */
00147   img->type=IMG_TYPE_IMAGE;
00148   strncpy(img->studyNr, dsr.hist.patient_id, MAX_STUDYNR_LEN);
00149   strcpy(img->patientName, dsr.hist.patient_id);
00150   img->sizex=dsr.dime.pixdim[1];
00151   img->sizey=dsr.dime.pixdim[2];
00152   img->sizez=dsr.dime.pixdim[3];
00153   /*if(dsr.dime.funused2>1.E-5) img->zoom=dsr.dime.funused2;*/
00154   if(dsr.dime.funused3>1.E-5) img->isotopeHalflife=dsr.dime.funused3;
00155   for(pi=0; pi<dimz; pi++) img->planeNumber[pi]=pi+1;
00156   if(dsr.little) img->_fileFormat=IMG_ANA_L; else img->_fileFormat=IMG_ANA;
00157   /* Decay correction */
00158   if(strstr(dsr.hist.descrip, "Decay corrected.")!=NULL)
00159     img->decayCorrected=1;
00160   else if(strstr(dsr.hist.descrip, "No decay correction.")!=NULL)
00161     img->decayCorrected=0;
00162   else
00163     img->decayCorrected=1;
00164 
00165   /* Allocate memory for one image frame */
00166   fdata=malloc(pxlNr*sizeof(float));
00167   if(fdata==NULL) {fclose(fp); imgSetStatus(img, STATUS_NOMEMORY); return(12);}
00168 
00169   /* Read one image frame at a time */
00170   for(fi=0; fi<dimt; fi++) {
00171     fptr=fdata;
00172     ret=anaReadImagedata(fp, &dsr, fi+1, fptr);
00173     if(ret) {free(fdata); fclose(fp); imgSetStatus(img, STATUS_NOIMGDATA); return(7);}
00174     /* Copy pixel values to IMG */
00175     fptr=fdata;
00176     if(anaFlipping()==0) { /* no flipping in z-direction */
00177       for(pi=0; pi<img->dimz; pi++)
00178         for(yi=dimy-1; yi>=0; yi--)
00179           for(xi=dimx-1; xi>=0; xi--)
00180             img->m[pi][yi][xi][fi]=*fptr++;
00181     } else {
00182       for(pi=dimz-1; pi>=0; pi--)
00183         for(yi=dimy-1; yi>=0; yi--)
00184           for(xi=dimx-1; xi>=0; xi--)
00185             img->m[pi][yi][xi][fi]=*fptr++;
00186     }
00187   } /* next frame */
00188   free(fdata);
00189   fclose(fp);
00190   
00191   /* Try to read frame time information from SIF file */
00192   /* Make filename from database or image data file */
00193   strcpy(siffile, dbname); strcat(siffile, ".sif");
00194   if(access(siffile, 0) == -1) {
00195     strcpy(siffile, datfile); strcat(siffile, ".sif");
00196   }
00197   /* Check if SIF file is found */
00198   if(IMG_TEST) printf("reading SIF file %s\n", siffile);
00199   if(access(siffile, 0) == -1) {
00200     if(IMG_TEST) printf(" No SIF file; therefore unknown frame times.\n");
00201     return(0);
00202   }
00203   /* If found, then read it */
00204   sifInit(&sif); ret=sifRead(siffile, &sif);
00205   if(ret) {imgSetStatus(img, STATUS_NOSIFDATA); return(21);}
00206   /* Copy scan start time */
00207   if(sif.scantime>0) {
00208     st=localtime(&sif.scantime);
00209     if(st!=NULL) strftime(buf, 128, "%Y-%m-%d %H:%M:%S", st);
00210     else strcpy(buf, "1900-01-01 00:00:00");
00211     img->scanStart=sif.scantime;
00212   }
00213   /* Copy frame times */
00214   if(sif.frameNr!=img->dimt) {
00215     imgSetStatus(img, STATUS_WRONGSIFDATA); sifEmpty(&sif); return(22);}
00216   for(fi=0; fi<sif.frameNr; fi++) {
00217     img->start[fi]=sif.x1[fi]; img->end[fi]=sif.x2[fi];
00218     img->mid[fi]=0.5*(img->start[fi]+img->end[fi]);
00219   }
00220   /* Copy prompts and randoms */
00221   for(fi=0; fi<sif.frameNr; fi++) {
00222     img->prompts[fi]=sif.prompts[fi];
00223     img->randoms[fi]=sif.randoms[fi];
00224   }
00225   
00226   /* Set isotopeHalflife, if isotope is found */
00227   if(strlen(sif.isotope_name)>1) {
00228     img->isotopeHalflife=60.0*hlFromIsotope(sif.isotope_name);
00229   }
00230   sifEmpty(&sif);
00231   
00232   return(0);
00233 }
00234 /*****************************************************************************/
00235 
00236 /*****************************************************************************/
00253 int imgWriteAnalyze(const char *dbname, IMG *img) {
00254   FILE *fp;
00255   int ret, fi, pi, xi, yi, little;
00256   float g;
00257   ANALYZE_DSR dsr;
00258   char datfile[FILENAME_MAX], hdrfile[FILENAME_MAX];
00259   const char *cptr;
00260   int pxlNr=0;
00261   struct tm *st;
00262   short int *sdata, *sptr, smin, smax;
00263 
00264 
00265   if(IMG_TEST) printf("imgWriteAnalyze(%s, *img)\n", dbname);
00266 
00267   /* Check the arguments */
00268   imgSetStatus(img, STATUS_OK);
00269   if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) {
00270     imgSetStatus(img, STATUS_FAULT); return(2);}
00271   if(dbname==NULL || !dbname[0]) {imgSetStatus(img, STATUS_FAULT); return(1);}
00272   
00273   /* Make the image and header filenames */  
00274   strcpy(datfile, dbname); strcat(datfile, ".img");
00275   strcpy(hdrfile, dbname); strcat(hdrfile, ".hdr");
00276 
00277 
00278   /*
00279    *  Fill Analyze header
00280    */
00281   if(img->_fileFormat==IMG_ANA_L) dsr.little=1; else dsr.little=0;
00282   /* Header key */
00283   memset(&dsr.hk, 0, sizeof(ANALYZE_HEADER_KEY));
00284   memset(&dsr.dime, 0, sizeof(ANALYZE_HEADER_IMGDIM));
00285   memset(&dsr.hist, 0, sizeof(ANALYZE_HEADER_HISTORY));
00286   dsr.hk.sizeof_hdr=348;
00287   strcpy(dsr.hk.data_type, "");
00288   cptr=strrchr(dbname, '/'); if(cptr==NULL) cptr=strrchr(dbname, '\\');
00289   if(cptr!=NULL) cptr++; if(cptr==NULL) cptr=dbname;
00290   strncpy(dsr.hk.db_name, cptr, 17);
00291   dsr.hk.extents=16384;
00292   dsr.hk.regular='r';
00293   /* Image dimension */
00294   dsr.dime.dim[0]=4;
00295   dsr.dime.dim[1]=img->dimx;
00296   dsr.dime.dim[2]=img->dimy;
00297   dsr.dime.dim[3]=img->dimz;
00298   dsr.dime.dim[4]=img->dimt;
00299   dsr.dime.datatype=ANALYZE_DT_SIGNED_SHORT;
00300   dsr.dime.bitpix=16;
00301   dsr.dime.pixdim[0]=0.0;
00302   dsr.dime.pixdim[1]=img->sizex;
00303   dsr.dime.pixdim[2]=img->sizey;
00304   dsr.dime.pixdim[3]=img->sizez;
00305   dsr.dime.pixdim[4]=0.0;
00306   dsr.dime.funused1=0.0; /* Scale factor is set later */
00307   /* dsr.dime.funused2=img->zoom; */ /* Reconstruction zoom */
00308   dsr.dime.funused3=img->isotopeHalflife;
00309   /* Data history */
00310   if(img->decayCorrected==1) strcpy(dsr.hist.descrip, "Decay corrected.");
00311   else strcpy(dsr.hist.descrip, "No decay correction.");
00312   strncpy(dsr.hist.scannum, img->studyNr, 10);
00313   st=localtime(&img->scanStart);
00314   if(st!=NULL) {
00315     strftime(dsr.hist.exp_date, 10, "%Y-%m-%d", st);
00316     strftime(dsr.hist.exp_time, 10, "%H:%M:%S", st);
00317   } else {
00318     strncpy(dsr.hist.exp_date, "1900-01-01", 10);
00319     strncpy(dsr.hist.exp_time, "00:00:00", 10);
00320   }
00321 
00322   /*
00323    *  Scale data to short int range
00324    *  Determine and set scale factor and cal_min & cal_max
00325    */
00326   if(IMG_TEST) printf("scaling data to short ints\n");
00327   ret=imgMinMax(img, &dsr.dime.cal_min, &dsr.dime.cal_max);
00328   if(ret) {imgSetStatus(img, STATUS_FAULT); return(3);}
00329   if(fabs(dsr.dime.cal_min)>fabs(dsr.dime.cal_max)) g=fabs(dsr.dime.cal_min);
00330   else g=fabs(dsr.dime.cal_max);
00331   if(g<1E-20) g=1.0; else g=32767./g; dsr.dime.funused1=1.0/g;
00332   if(IMG_TEST) printf("min=%g max=%g scale_factor=%g\n",
00333     dsr.dime.cal_min, dsr.dime.cal_max, dsr.dime.funused1);
00334 
00335   /* Allocate memory for short int array */
00336   pxlNr=(img->dimx)*(img->dimy)*(img->dimz);
00337   sdata=malloc(pxlNr*sizeof(short int));
00338   if(sdata==NULL) {
00339     imgSetStatus(img, STATUS_NOMEMORY);
00340     return 12;
00341   }
00342 
00343   /* Open image data file for write */
00344   if((fp=fopen(datfile, "wb")) == NULL) {
00345     imgSetStatus(img, STATUS_CANTWRITEIMGFILE);
00346     free(sdata);
00347     return 14;
00348   }
00349 
00350   /* Copy and write image matrix data to short int array */
00351   /* Data is written one frame at a time */
00352   smin=smax=temp_roundf(g*img->m[0][0][0][0]);
00353   for(fi=0; fi<img->dimt; fi++) {
00354     sptr=sdata;
00355     if(anaFlipping()==0) {
00356       for(pi=0; pi<img->dimz; pi++)
00357         for(yi=img->dimy-1; yi>=0; yi--)
00358           for(xi=img->dimx-1; xi>=0; xi--) {
00359             *sptr=temp_roundf(g*img->m[pi][yi][xi][fi]);
00360             if(*sptr>smax) smax=*sptr; else if(*sptr<smin) smin=*sptr;
00361             sptr++;
00362           }
00363     } else {
00364       for(pi=img->dimz-1; pi>=0; pi--)
00365         for(yi=img->dimy-1; yi>=0; yi--)
00366           for(xi=img->dimx-1; xi>=0; xi--) {
00367             *sptr=temp_roundf(g*img->m[pi][yi][xi][fi]);
00368             if(*sptr>smax) smax=*sptr; else if(*sptr<smin) smin=*sptr;
00369             sptr++;
00370           }
00371     }
00372     /* Change byte order if necessary */
00373     little=little_endian();
00374     if(little!=dsr.little)
00375       swabip(sdata, pxlNr*sizeof(short int));
00376     /* Write image data */
00377     if(fwrite(sdata, 2, pxlNr, fp) != pxlNr) {
00378       imgSetStatus(img, STATUS_CANTWRITEIMGFILE);
00379       free(sdata); fclose(fp);
00380       return 15;
00381     }
00382   }
00383   /* Done writing */
00384   fclose(fp);
00385   free(sdata);
00386 
00387   if(IMG_TEST) printf("smin=%d smax=%d\n", smin, smax);
00388 
00389   /* Set header glmin & glmax */
00390   dsr.dime.glmin=smin; dsr.dime.glmax=smax;
00391   
00392   /* Write Analyze header */
00393   ret=anaWriteHeader(hdrfile, &dsr);
00394   if(ret) {
00395     imgSetStatus(img, STATUS_CANTWRITEHEADERFILE);
00396     return 21;
00397   }
00398   imgSetStatus(img, STATUS_OK);
00399   return 0;
00400 }
00401 /*****************************************************************************/
00402 
00403 /*****************************************************************************/
00414 int imgReadAnalyzeHeader(const char *dbname, IMG *img) {
00415   char hdrfile[FILENAME_MAX], siffile[FILENAME_MAX];
00416   ANALYZE_DSR ana_header;
00417   SIF sif;
00418   double f;
00419   int ret;
00420 
00421   if(IMG_TEST) printf("\nimgReadAnalyzeHeader(%s, *img)\n", dbname);
00422   
00423   /* Check the input */
00424   if(img==NULL) return STATUS_FAULT;
00425   if(img->status!=IMG_STATUS_INITIALIZED) return STATUS_FAULT;
00426   imgSetStatus(img, STATUS_FAULT);
00427   if(dbname==NULL) return STATUS_FAULT;
00428 
00429   /* Determine the names of hdr and sif files */
00430   ret=anaDatabaseExists(dbname, hdrfile, NULL, siffile);
00431   if(ret==0) return STATUS_NOFILE;
00432   
00433   /* Read Analyze header file */
00434   ret=anaReadHeader(hdrfile, &ana_header);
00435   if(ret!=0) {
00436     if(IMG_TEST>1) printf("anaReadHeader() return value := %d\n", ret);
00437     if(ret==1) return STATUS_FAULT;
00438     else if(ret==2) return STATUS_NOHEADERFILE;
00439     else return STATUS_UNSUPPORTED;
00440     return(STATUS_FAULT);
00441   }
00442   /* and set IMG contents */
00443   ret=imgGetAnalyzeHeader(img, &ana_header);
00444   if(ret!=0) {
00445     imgSetStatus(img, ret);
00446     return(ret);
00447   }
00448 
00449   /* If SIF does not exist, then that's it */
00450   if(!siffile[0]) {
00451     imgSetStatus(img, STATUS_OK);
00452     return STATUS_OK;
00453   }
00454 
00455   /* SIF is available, so read that too */
00456   sifInit(&sif); ret=0;
00457   if(sifRead(siffile, &sif)!=0) return STATUS_OK;
00458   /* Copy scan time */
00459   img->scanStart=sif.scantime;
00460   /* Study number, if not yet defined */
00461   if(!img->studyNr[0] && strlen(sif.studynr)>1 )
00462     strncpy(img->studyNr, sif.studynr, MAX_STUDYNR_LEN);
00463   /* Isotope half-life, if not yet defined */
00464   f=hlFromIsotope(sif.isotope_name);
00465   if(img->isotopeHalflife<=0.0 && f>0.0) img->isotopeHalflife=60.0*f;
00466   sifEmpty(&sif);
00467 
00468   return STATUS_OK;
00469 }
00470 /*****************************************************************************/
00471 
00472 /*****************************************************************************/
00481 int imgGetAnalyzeHeader(IMG *img, ANALYZE_DSR *h) {
00482   int dimNr, dimx, dimy, dimz=1, dimt=1, pxlNr=0;
00483 
00484   if(IMG_TEST) printf("\nimgGetAnalyzeHeader(*img, *dsr)\n");
00485   
00486   /* Check the input */
00487   
00488   if(img==NULL) return STATUS_FAULT;
00489   if(img->status!=IMG_STATUS_INITIALIZED && img->status!=IMG_STATUS_OCCUPIED)
00490     return STATUS_FAULT;
00491   imgSetStatus(img, STATUS_FAULT);
00492   if(h==NULL) return STATUS_FAULT;
00493     
00494   imgSetStatus(img, STATUS_INVALIDHEADER);
00495 
00496   /* Get the image dimensions from header */
00497   dimNr=h->dime.dim[0];
00498   if(dimNr<2) return STATUS_INVALIDHEADER;
00499   dimx=h->dime.dim[1]; dimy=h->dime.dim[2];
00500   if(dimNr>2) {dimz=h->dime.dim[3]; if(dimNr>3) dimt=h->dime.dim[4];}
00501   pxlNr=dimx*dimy*dimz;
00502   if(pxlNr<1) return STATUS_INVALIDHEADER;
00503   img->dimx=dimx; img->dimy=dimy; img->dimz=dimz; img->dimt=dimt;
00504 
00505   /* Copy information from Analyze header */
00506   img->type=IMG_TYPE_IMAGE;
00507   strncpy(img->studyNr, h->hist.patient_id, MAX_STUDYNR_LEN);
00508   strcpy(img->patientName, h->hist.patient_id);
00509   img->sizex=h->dime.pixdim[1];
00510   img->sizey=h->dime.pixdim[2];
00511   img->sizez=h->dime.pixdim[3];
00512   /*if(h->dime.funused2>1.E-5) img->zoom=h->dime.funused2;*/
00513   if(h->dime.funused3>1.E-5) img->isotopeHalflife=h->dime.funused3;
00514   if(h->little) img->_fileFormat=IMG_ANA_L; else img->_fileFormat=IMG_ANA;
00515   
00516   /* Decay correction */
00517   if(strstr(h->hist.descrip, "Decay corrected.")!=NULL)
00518     img->decayCorrected=1;
00519   else if(strstr(h->hist.descrip, "No decay correction.")!=NULL)
00520     img->decayCorrected=0;
00521   else
00522     img->decayCorrected=1;
00523 
00524   imgSetStatus(img, STATUS_OK);
00525   return STATUS_OK;
00526 }
00527 /*****************************************************************************/
00528 
00529 /*****************************************************************************/
00542 int imgSetAnalyzeHeader(IMG *img, const char *dbname, ANALYZE_DSR *dsr, float fmin, float fmax) {
00543   struct tm *st;
00544   char *cptr;
00545   float g;
00546 
00547   if(IMG_TEST) printf("\nimgSetAnalyzeHeader(*img, *dsr)\n");
00548   
00549   /* Check the input */
00550   if(img==NULL) return STATUS_FAULT;
00551   if(img->status!=IMG_STATUS_INITIALIZED && img->status!=IMG_STATUS_OCCUPIED)
00552     return STATUS_FAULT;
00553   imgSetStatus(img, STATUS_FAULT);
00554   if(dsr==NULL) return STATUS_FAULT;
00555 
00556   /* Byte order */
00557   if(img->_fileFormat==IMG_ANA_L) dsr->little=1; else dsr->little=0;
00558   /* Header key */
00559   memset(&dsr->hk, 0, sizeof(ANALYZE_HEADER_KEY));
00560   memset(&dsr->dime, 0, sizeof(ANALYZE_HEADER_IMGDIM));
00561   memset(&dsr->hist, 0, sizeof(ANALYZE_HEADER_HISTORY));
00562   dsr->hk.sizeof_hdr=348;
00563   strcpy(dsr->hk.data_type, "");
00564   cptr=strrchr(dbname, '/'); if(cptr==NULL) cptr=strrchr(dbname, '\\');
00565   if(cptr!=NULL) cptr++; if(cptr==NULL) cptr=(char*)dbname;
00566   strncpy(dsr->hk.db_name, cptr, 17);
00567   dsr->hk.extents=16384;
00568   dsr->hk.regular='r';
00569   /* Image dimension */
00570   dsr->dime.dim[0]=4;
00571   dsr->dime.dim[1]=img->dimx;
00572   dsr->dime.dim[2]=img->dimy;
00573   dsr->dime.dim[3]=img->dimz;
00574   dsr->dime.dim[4]=img->dimt;
00575   dsr->dime.datatype=ANALYZE_DT_SIGNED_SHORT;
00576   dsr->dime.bitpix=16;
00577   dsr->dime.pixdim[0]=0.0;
00578   dsr->dime.pixdim[1]=img->sizex;
00579   dsr->dime.pixdim[2]=img->sizey;
00580   dsr->dime.pixdim[3]=img->sizez;
00581   dsr->dime.pixdim[4]=0.0;
00582   dsr->dime.funused1=0.0; /* Scale factor is set later */
00583   /* dsr.dime.funused2=img->zoom; */ /* Reconstruction zoom */
00584   dsr->dime.funused3=img->isotopeHalflife;
00585   /* Data history */
00586   if(img->decayCorrected==1) strcpy(dsr->hist.descrip, "Decay corrected.");
00587   else strcpy(dsr->hist.descrip, "No decay correction.");
00588   strncpy(dsr->hist.scannum, img->studyNr, 10);
00589   st=localtime(&img->scanStart);
00590   if(st!=NULL) {
00591     strftime(dsr->hist.exp_date, 10, "%Y-%m-%d", st);
00592     strftime(dsr->hist.exp_time, 10, "%H:%M:%S", st);
00593   } else {
00594     strcpy(dsr->hist.exp_date, "1900-01-01");
00595     strcpy(dsr->hist.exp_time, "00:00:00");
00596   }
00597   /* Determine and set scale factor and cal_min & cal_max */
00598   if(fmin<fmax) {
00599     dsr->dime.cal_min=fmin; dsr->dime.cal_max=fmax;
00600   } else { /* not given in function call, try to find those here */
00601     if(img->status==IMG_STATUS_OCCUPIED &&
00602        imgMinMax(img, &dsr->dime.cal_min, &dsr->dime.cal_max)==0) {}
00603     else return STATUS_FAULT;
00604   }
00605   if(fabs(dsr->dime.cal_min) > fabs(dsr->dime.cal_max)) g = fabs(dsr->dime.cal_min);
00606   else g = fabs(dsr->dime.cal_max);
00607   /* if(fabs(dsr->dime.cal_min)>fabs(dsr->dime.cal_max)) g=fabs(dsr->dime.cal_min); */
00608   /* else g=fabs(dsr->dime.cal_max); */
00609   if(g<1E-20) g=1.0; else g=32767./g; dsr->dime.funused1=1.0/g;
00610   /* Set header glmin & glmax */
00611   dsr->dime.glmin=temp_roundf(fmin*g); dsr->dime.glmax=temp_roundf(fmax*g);
00612   /* printf("glmin=%d\n", dsr->dime.glmin); */
00613   /* printf("glmax=%d\n", dsr->dime.glmax); */
00614 
00615   imgSetStatus(img, STATUS_OK);
00616   return STATUS_OK;
00617 }
00618 /*****************************************************************************/
00619 
00620 /*****************************************************************************/
00629 int imgReadAnalyzeFirstFrame(const char *fname, IMG *img) {
00630   int ret=0;
00631 
00632   if(IMG_TEST) printf("\nimgReadAnalyzeFirstFrame(%s, *img)\n", fname);
00633   /* Check the input */
00634   if(img==NULL) return STATUS_FAULT;
00635   if(img->status!=IMG_STATUS_INITIALIZED) return STATUS_FAULT;
00636   imgSetStatus(img, STATUS_FAULT);
00637   if(fname==NULL) return STATUS_FAULT;
00638 
00639   /* Read header information from file */
00640   ret=imgReadAnalyzeHeader(fname, img); if(ret) return(ret);
00641   if(IMG_TEST>3) imgInfo(img);
00642 
00643   /* Allocate memory for one frame */
00644   img->dimt=1;
00645   ret=imgAllocate(img, img->dimz, img->dimy, img->dimx, img->dimt);
00646   if(ret) return STATUS_NOMEMORY;
00647 
00648   /* Read the first frame */
00649   ret=imgReadAnalyzeFrame(fname, 1, img, 0); if(ret) return(ret); 
00650 
00651   /* All went well */
00652   imgSetStatus(img, STATUS_OK);
00653   return STATUS_OK;
00654 }
00655 /*****************************************************************************/
00656 
00657 /*****************************************************************************/
00674 int imgReadAnalyzeFrame(const char *fname, int frame_to_read, IMG *img, int frame_index) {
00675   FILE *fp;
00676   int ret, pi, xi, yi;
00677   float *fdata=NULL, *fptr;
00678   ANALYZE_DSR dsr;
00679   char datfile[FILENAME_MAX], hdrfile[FILENAME_MAX], siffile[FILENAME_MAX];
00680   SIF sif;
00681 
00682 
00683   if(IMG_TEST) printf("\nimgReadAnalyzeFrame(%s, %d, *img, %d)\n",
00684     fname, frame_to_read, frame_index);
00685     
00686   /* Check the input */
00687   if(img==NULL) return STATUS_FAULT;
00688   if(img->status!=IMG_STATUS_OCCUPIED) return STATUS_FAULT;
00689   if(fname==NULL) return STATUS_FAULT;
00690   if(frame_index<0 || frame_index>img->dimt-1) return STATUS_FAULT;
00691   if(frame_to_read<1) return STATUS_FAULT;
00692   imgSetStatus(img, STATUS_FAULT);
00693   
00694   /* Determine the names of hdr, data and sif files */
00695   ret=anaDatabaseExists(fname, hdrfile, datfile, siffile);
00696   if(ret==0) return STATUS_NOFILE;
00697 
00698   /* Read Analyze header file */
00699   ret=anaReadHeader(hdrfile, &dsr);
00700   if(ret!=0) {
00701     if(ret==1) return STATUS_FAULT;
00702     else if(ret==2) return STATUS_NOHEADERFILE;
00703     else return STATUS_UNSUPPORTED;
00704     return(STATUS_FAULT);
00705   }
00706 
00707   /* Open image datafile */
00708   if(IMG_TEST>2) fprintf(stdout, "reading image data %s\n", datfile);
00709   if((fp=fopen(datfile, "rb")) == NULL) return STATUS_NOIMGDATA;
00710 
00711   /* Allocate memory for one image frame */
00712   fdata=malloc(img->dimx*img->dimy*img->dimz*sizeof(float));
00713   if(fdata==NULL) {fclose(fp); return STATUS_NOMEMORY;}
00714 
00715   /* Read the required image frame */
00716   fptr=fdata;
00717   ret=anaReadImagedata(fp, &dsr, frame_to_read, fptr);
00718   fclose(fp);
00719   if(ret==3) {free(fdata); return STATUS_NOMATRIX;} /* no more frames */
00720   if(ret!=0) {free(fdata); return STATUS_UNSUPPORTED;}
00721 
00722   /* Copy pixel values to IMG */
00723   fptr=fdata;
00724   if(anaFlipping()==0) { /* no flipping in z-direction */
00725     for(pi=0; pi<img->dimz; pi++)
00726       for(yi=img->dimy-1; yi>=0; yi--)
00727         for(xi=img->dimx-1; xi>=0; xi--)
00728           img->m[pi][yi][xi][frame_index]=*fptr++;
00729   } else {
00730     for(pi=img->dimz-1; pi>=0; pi--)
00731       for(yi=img->dimy-1; yi>=0; yi--)
00732         for(xi=img->dimx-1; xi>=0; xi--)
00733           img->m[pi][yi][xi][frame_index]=*fptr++;
00734   }
00735   free(fdata);
00736 
00737   imgSetStatus(img, STATUS_OK); /* If the rest is failed, no problem */
00738 
00739   /* 
00740    *  Try to read frame time information from SIF file
00741    */
00742   sifInit(&sif);
00743   if(sifRead(siffile, &sif)!=0) return STATUS_OK;
00744   /* Frame information */
00745   if(sif.frameNr>=frame_to_read) {
00746     img->start[frame_index]=sif.x1[frame_to_read-1];
00747     img->end[frame_index]=sif.x2[frame_to_read-1];
00748     img->mid[frame_index]=0.5*(img->start[frame_index]+img->end[frame_index]);
00749     img->prompts[frame_index]=sif.prompts[frame_to_read-1];
00750     img->randoms[frame_index]=sif.randoms[frame_to_read-1];
00751   }
00752   sifEmpty(&sif);
00753 
00754   return STATUS_OK;
00755 }
00756 /*****************************************************************************/
00757 
00758 /*****************************************************************************/
00781 int imgWriteAnalyzeFrame(const char *dbname, int frame_to_write, IMG *img, int frame_index, float fmin, float fmax) {
00782   IMG test_img;
00783   int ret=0, pxlNr, zi, xi, yi, little;
00784   FILE *fp;
00785   short int *sdata=NULL, *sptr;
00786   char datfile[FILENAME_MAX], hdrfile[FILENAME_MAX], siffile[FILENAME_MAX];
00787   ANALYZE_DSR dsr;
00788   float scale_factor=1.0;
00789 
00790 
00791   if(IMG_TEST) printf("\nimgWriteAnalyzeFrame(%s, %d, *img, %d, %g, %g)\n",
00792     dbname, frame_to_write, frame_index, fmin, fmax);
00793 
00794   /*
00795    *  Check the input 
00796    */
00797   if(dbname==NULL) return STATUS_FAULT;
00798   if(img==NULL) return STATUS_FAULT;
00799   if(img->status!=IMG_STATUS_OCCUPIED) return STATUS_FAULT;
00800   if(frame_to_write<0) return STATUS_FAULT;
00801   if(frame_index<0 || frame_index>=img->dimt) return STATUS_FAULT;
00802   if(img->_fileFormat!=IMG_ANA_L && img->_fileFormat!=IMG_ANA) return STATUS_FAULT;
00803 
00804   /*
00805    *  If database does not exist, then create it with new header,
00806    *  and if it does exist, then read and check header information.
00807    *  Create or edit header to contain correct frame nr.
00808    *  Determine the global scaling factor.   
00809    */
00810   imgInit(&test_img);
00811   if(anaDatabaseExists(dbname, hdrfile, datfile, siffile)==0) { /* not existing */
00812 
00813     /* Create database filenames */
00814     sprintf(hdrfile, "%s.hdr", dbname);
00815     sprintf(datfile, "%s.img", dbname);
00816     sprintf(siffile, "%s.sif", dbname);
00817 
00818     /* Set main header */
00819     imgSetAnalyzeHeader(img, dbname, &dsr, fmin, fmax);
00820     if(frame_to_write==0) frame_to_write=1;
00821     dsr.dime.dim[4]=frame_to_write;
00822     scale_factor=dsr.dime.funused1;
00823     if(fabs(scale_factor)>1.0E-20) scale_factor=1.0/scale_factor;
00824 
00825     /* Write Analyze header */
00826     ret=anaWriteHeader(hdrfile, &dsr); if(ret && IMG_TEST) printf("anaWriteHeader() := %d\n", ret);
00827     if(ret) return STATUS_CANTWRITEHEADERFILE;
00828 
00829     /* Remove datafile if necessary */
00830     if(access(datfile, 0) != -1) remove(datfile);
00831 
00832   } else { /* database does exist */
00833   
00834     /* Read header information for checking */
00835     ret=imgReadAnalyzeHeader(dbname, &test_img); if(ret!=0) return ret;
00836     /* Check that file format is the same */
00837     if(img->_fileFormat!=test_img._fileFormat || img->type!=test_img.type)
00838       return STATUS_WRONGFILETYPE;
00839     /* Check that matrix sizes are the same */
00840     if(img->dimz!=test_img.dimz || img->dimx!=test_img.dimx ||
00841        img->dimy!=test_img.dimy)
00842       return STATUS_VARMATSIZE;
00843     imgEmpty(&test_img);
00844 
00845     /* Read the header, set new frame number, and write it back */
00846     /* Get also the scale factor */
00847     if((ret=anaReadHeader(hdrfile, &dsr))!=0) return STATUS_NOMAINHEADER;
00848     scale_factor=1.0/dsr.dime.funused1;
00849     if(frame_to_write==0) frame_to_write=dsr.dime.dim[4]+1;
00850     if(dsr.dime.dim[4]<frame_to_write) {
00851       if(dsr.dime.dim[4]+1<frame_to_write) return STATUS_MISSINGMATRIX;
00852       dsr.dime.dim[4]=frame_to_write;
00853     }
00854     if((ret=anaWriteHeader(hdrfile, &dsr))!=0) return STATUS_NOWRITEPERM;
00855   }
00856   if(IMG_TEST>2) {
00857     printf("frame_to_write := %d\n", frame_to_write);
00858     printf("hdrfile := %s\n", hdrfile);
00859     printf("datfile := %s\n", datfile);
00860     printf("siffile := %s\n", siffile);
00861   }
00862 
00863   /* Allocate memory for matrix short int data (one plane) */
00864   pxlNr=img->dimx*img->dimy;
00865   sdata=(short int*)malloc(pxlNr*sizeof(short int));
00866   if(sdata==NULL) return STATUS_NOMEMORY;
00867 
00868   /* Open datafile, not removing possible old contents */
00869   if(frame_to_write==1) fp=fopen(datfile, "wb"); else fp=fopen(datfile, "r+b");
00870   if(fp==NULL) {free(sdata); return STATUS_CANTWRITEIMGFILE;}
00871   /* Move file pointer to the place of current frame */
00872   if(fseek(fp, (frame_to_write-1)*pxlNr*img->dimz*sizeof(short int), SEEK_SET)!=0) {
00873     free(sdata); fclose(fp); return STATUS_MISSINGMATRIX;}
00874   little=little_endian();
00875   /* Copy, scale and write data plane-by-plane */
00876   if(anaFlipping()==0) {
00877     for(zi=0; zi<img->dimz; zi++) {
00878       sptr=sdata; /*printf("plane := %d\n  scale_factor := %g\n", zi+1, scale_factor);*/
00879       for(yi=img->dimy-1; yi>=0; yi--) for(xi=img->dimx-1; xi>=0; xi--) {
00880         *sptr=temp_roundf(scale_factor*img->m[zi][yi][xi][frame_index]); sptr++;
00881       }
00882       /* Change byte order if necessary */
00883       sptr=sdata; if(little!=dsr.little) swabip(sptr, pxlNr*sizeof(short int));
00884       /* Write image data */
00885       sptr=sdata;
00886       if(fwrite(sptr, sizeof(short int), pxlNr, fp) != pxlNr) {
00887         free(sdata); fclose(fp); return STATUS_CANTWRITEIMGFILE;
00888       }
00889     }
00890   } else {
00891     for(zi=img->dimz-1; zi>=0; zi--) {
00892       sptr=sdata;
00893       for(yi=img->dimy-1; yi>=0; yi--) for(xi=img->dimx-1; xi>=0; xi--) {
00894         *sptr=temp_roundf(scale_factor*img->m[zi][yi][xi][frame_index]); sptr++;
00895       }
00896       /* Change byte order if necessary */
00897       sptr=sdata; if(little!=dsr.little) swabip(sptr, pxlNr*sizeof(short int));
00898       /* Write image data */
00899       sptr=sdata;
00900       if(fwrite(sptr, sizeof(short int), pxlNr, fp) != pxlNr) {
00901         free(sdata); fclose(fp); return STATUS_CANTWRITEIMGFILE;
00902       }
00903     }
00904   }
00905   free(sdata);
00906   fclose(fp);
00907 
00908   return STATUS_OK;
00909 }
00910 /*****************************************************************************/
00911 
00912 /*****************************************************************************/