img_e63.c
Go to the documentation of this file.
00001 /******************************************************************************
00002 
00003   Copyright (c) 2007,2008 Turku PET Centre
00004 
00005   Library:     img_e63.c
00006   Description: ECAT 6.3 I/O routines for IMG data.
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     If valid study number is found in user_process_code, then take it from there.
00029     Patient_id and study_description are read and written.
00030     Prompts and randoms (delayed) are read and written.
00031   2007-03-25 VO
00032     Added functions imgGetEcat63MHeader(), imgSetEcat63MHeader(),
00033     imgEcat63Supported(), imgGetEcat63Fileformat(), imgReadEcat63Header(),
00034     imgReadEcat63FirstFrame(), imgReadEcat63Frame(), imgWriteEcat63Frame(),
00035     and imgSetEcat63SHeader().
00036   2007-09-10 VO
00037     Return value of localtime() is checked.
00038   2008-11-06 VO
00039     study_description is copied with strncpy(), not strcpy().
00040 
00041 
00042 
00043 ******************************************************************************/
00044 #include <stdio.h>
00045 #include <stdlib.h>
00046 #include <unistd.h>
00047 #include <math.h>
00048 #include <string.h>
00049 #include <time.h>
00050 /*****************************************************************************/
00051 #include "petc99.h"
00052 #include "swap.h"
00053 #include "halflife.h"
00054 /*****************************************************************************/
00055 #include "include/img.h"
00056 #include "include/ecat63.h"
00057 #include "include/ecat7.h"
00058 #include "include/imgmax.h"
00059 #include "include/imgdecay.h"
00060 #include "include/sif.h"
00061 #include "include/imgfile.h"
00062 /*****************************************************************************/
00063 
00064 /*****************************************************************************/
00077 int ecat63ReadAllToImg(const char *fname, IMG *img) {
00078   int                 i, j, m, ret, blkNr=0, dim_x=0, dim_y=0, pxlNr=0;
00079   int                 frameNr, planeNr, del_nr=0;
00080   int                 frame, plane, prev_frame, prev_plane, seqplane;
00081   FILE               *fp;
00082   ECAT63_mainheader   main_header;
00083   MATRIXLIST          mlist;
00084   MatDir              tmpMatdir;
00085   Matval              matval;
00086   ECAT63_imageheader  image_header;
00087   ECAT63_scanheader   scan_header;
00088   ECAT63_attnheader   attn_header;
00089   ECAT63_normheader   norm_header;
00090   float               scale=1.0;
00091   short int          *sptr;
00092   char               *mdata=NULL, *mptr;
00093   int                *iptr;
00094   struct tm           scanStart;
00095 
00096 
00097   if(IMG_TEST) printf("ecat63ReadAllToImg(%s, *img)\n", fname);
00098   /* Check the arguments */
00099   if(fname==NULL || img==NULL || img->status!=IMG_STATUS_INITIALIZED) {
00100     if(img) imgSetStatus(img, STATUS_FAULT);
00101     return 1;
00102   }
00103 
00104   /* Open input CTI file for read */
00105   if((fp=fopen(fname, "rb")) == NULL) {
00106     imgSetStatus(img, STATUS_NOFILE);
00107     return 3;
00108   }
00109 
00110   /* Read main header */
00111   if((ret=ecat63ReadMainheader(fp, &main_header))) {
00112     imgSetStatus(img, STATUS_NOMAINHEADER);
00113     return 4;
00114   }
00115   if(IMG_TEST>5) ecat63PrintMainheader(&main_header, stdout);
00116 
00117   /* Read matrix list and nr */
00118   ecat63InitMatlist(&mlist);
00119   ret=ecat63ReadMatlist(fp, &mlist);
00120   if(ret) {
00121     imgSetStatus(img, STATUS_NOMATLIST);
00122     return 5;
00123   }
00124   if(mlist.matrixNr<=0) {
00125     strcpy(ecat63errmsg, "matrix list is empty");
00126     return 5;
00127   }
00128   /* Sort matrix list */
00129   for(i=0; i<mlist.matrixNr-1; i++) for(j=i+1; j<mlist.matrixNr; j++) {
00130     if(mlist.matdir[i].matnum>mlist.matdir[j].matnum) {
00131       tmpMatdir=mlist.matdir[i];
00132       mlist.matdir[i]=mlist.matdir[j]; mlist.matdir[j]=tmpMatdir;
00133     }
00134   }
00135   if(IMG_TEST>6) {
00136     printf("matrix list after sorting:\n");
00137     ecat63PrintMatlist(&mlist);
00138   }
00139 
00140   /* Trim extra frames */
00141   if(main_header.num_frames>0) {
00142     del_nr=ecat63DeleteLateFrames(&mlist, main_header.num_frames);
00143     if(IMG_TEST && del_nr>0)
00144       printf("  %d entries in matrix list will not be used.\n", del_nr);
00145   }
00146   /* Calculate the number of planes, frames and gates */
00147   /* Check that all planes have equal nr of frames (gates) */
00148   /* and check that frames (gates) are consequentally numbered */
00149   prev_plane=plane=-1; prev_frame=frame=-1; frameNr=planeNr=0; ret=0;
00150   for(m=0; m<mlist.matrixNr; m++) if(mlist.matdir[m].matstat==1) {
00151     /* get frame and plane */
00152     mat_numdoc(mlist.matdir[m].matnum, &matval);
00153     plane=matval.plane;
00154     if(main_header.num_frames>=main_header.num_gates) frame=matval.frame;
00155     else frame=matval.gate;
00156     if(IMG_TEST>2) printf("frame=%d plane=%d\n", frame, plane);
00157     /* did plane change? */
00158     if(plane!=prev_plane) {
00159       frameNr=1; planeNr++;
00160     } else {
00161       frameNr++;
00162       /* In each plane, frame (gate) numbers must be continuous */
00163       if(prev_frame>0 && frame!=prev_frame+1) {ret=1; break;}
00164     }
00165     prev_plane=plane; prev_frame=frame;
00166     /* Calculate and check the size of one data matrix */
00167     if(m==0) {
00168       blkNr=mlist.matdir[m].endblk-mlist.matdir[m].strtblk;
00169     } else if(blkNr!=mlist.matdir[m].endblk-mlist.matdir[m].strtblk) {
00170       ret=2; break;
00171     }
00172   } /* next matrix */
00173   if(IMG_TEST) printf("frameNr=%d planeNr=%d\n", frameNr, planeNr);
00174   if(ret==1 || (frameNr*planeNr != mlist.matrixNr-del_nr)) {
00175     strcpy(ecat63errmsg, "missing matrix");
00176     ecat63EmptyMatlist(&mlist); fclose(fp);
00177     return(6); /* this number is used in imgRead() */
00178   }
00179   if(ret==2) {
00180     strcpy(ecat63errmsg, "matrix sizes are different");
00181     ecat63EmptyMatlist(&mlist); fclose(fp); return(7);
00182   }
00183   /* Read x,y-dimensions from 1st matrix subheader */
00184   m=0; if(main_header.file_type==IMAGE_DATA) {
00185     ret=ecat63ReadImageheader(fp, mlist.matdir[m].strtblk, &image_header);
00186     dim_x=image_header.dimension_1; dim_y=image_header.dimension_2;
00187   } else if(main_header.file_type==RAW_DATA) {
00188     ret=ecat63ReadScanheader(fp, mlist.matdir[m].strtblk, &scan_header);
00189     dim_x=scan_header.dimension_1; dim_y=scan_header.dimension_2;
00190   } else if(main_header.file_type==ATTN_DATA) {
00191     ret=ecat63ReadAttnheader(fp, mlist.matdir[m].strtblk, &attn_header);
00192     dim_x=attn_header.dimension_1; dim_y=attn_header.dimension_2;
00193   } else if(main_header.file_type==NORM_DATA) {
00194     ret=ecat63ReadNormheader(fp, mlist.matdir[m].strtblk, &norm_header);
00195     dim_x=norm_header.dimension_1; dim_y=norm_header.dimension_2;
00196   }
00197   pxlNr=dim_x*dim_y;
00198   if(ret) {
00199     sprintf(ecat63errmsg, "cannot read matrix %u subheader", mlist.matdir[m].matnum);
00200     ecat63EmptyMatlist(&mlist); fclose(fp); return(8);
00201   }
00202 
00203   /* Allocate memory for ECAT data matrix */
00204   if(IMG_TEST>1) printf("allocating memory for %d blocks\n", blkNr);
00205   mdata=(char*)malloc(blkNr*MatBLKSIZE); if(mdata==NULL) {
00206     strcpy(ecat63errmsg, "out of memory");
00207     fclose(fp); ecat63EmptyMatlist(&mlist); return(8);
00208   }
00209   /* Allocate memory for whole img data */
00210   ret=imgAllocate(img, planeNr, dim_y, dim_x, frameNr);
00211   if(ret) {
00212     sprintf(ecat63errmsg, "out of memory (%d)", ret);
00213     fclose(fp); ecat63EmptyMatlist(&mlist); free(mdata); return(9);
00214   }
00215 
00216   /* Fill img info with ECAT main and sub header information */
00217   img->scanner=main_header.system_type;
00218   img->unit=main_header.calibration_units; /* this continues below */
00219   strncpy(img->radiopharmaceutical, main_header.radiopharmaceutical, 32);
00220   img->isotopeHalflife=main_header.isotope_halflife;
00221   memset(&scanStart, 0, sizeof(struct tm));
00222   scanStart.tm_year=main_header.scan_start_year-1900;
00223   scanStart.tm_mon=main_header.scan_start_month-1;
00224   scanStart.tm_mday=main_header.scan_start_day; scanStart.tm_yday=0;
00225   scanStart.tm_hour=main_header.scan_start_hour;
00226   scanStart.tm_min=main_header.scan_start_minute;
00227   scanStart.tm_sec=main_header.scan_start_second;
00228   scanStart.tm_isdst=-1;
00229   img->scanStart=mktime(&scanStart); if(img->scanStart==-1) img->scanStart=0;
00230   img->axialFOV=10.*main_header.axial_fov;
00231   img->transaxialFOV=10.*main_header.transaxial_fov;
00232   strncpy(img->studyNr, main_header.study_name, MAX_STUDYNR_LEN);
00233   strncpy(img->patientName, main_header.patient_name, 31);
00234   strncpy(img->patientID, main_header.patient_id, 15);
00235   img->sizez=10.*main_header.plane_separation;
00236   if(main_header.file_type==IMAGE_DATA) {
00237     img->type=IMG_TYPE_IMAGE;
00238     if(img->unit<1) img->unit=image_header.quant_units;
00239     img->zoom=image_header.recon_scale;
00240     if(image_header.decay_corr_fctr>1.0) img->decayCorrected=1;
00241     img->sizex=img->sizey=10.*image_header.pixel_size;
00242     if(img->sizez<10.*image_header.slice_width)
00243       img->sizez=10.*image_header.slice_width;
00244   } else if(main_header.file_type==RAW_DATA) {
00245     img->type=IMG_TYPE_RAW;
00246     img->decayCorrected=0;
00247   } else if(main_header.file_type==ATTN_DATA) {
00248     img->type=IMG_TYPE_RAW;
00249     img->decayCorrected=0;
00250   } else if(main_header.file_type==NORM_DATA) {
00251     img->type=IMG_TYPE_RAW;
00252     img->decayCorrected=0;
00253   }
00254   strncpy(img->studyDescription, main_header.study_description, 32);
00255   strncpy(img->userProcessCode, main_header.user_process_code, 10);
00256   img->userProcessCode[10]=(char)0;
00257   /* If valid study number is found in user_process_code, then take it */  
00258   if(!img->studyNr[0] && studynr_validity_check(img->userProcessCode))
00259     strcpy(img->studyNr, img->userProcessCode);
00260     
00261   /* Set _fileFormat */
00262   img->_fileFormat=IMG_E63;
00263 
00264   /* Read one ECAT matrix at a time and put them to img */
00265   prev_plane=plane=-1; seqplane=-1;
00266   for(m=0; m<mlist.matrixNr; m++) if(mlist.matdir[m].matstat==1) {
00267     /* get plane */
00268     mat_numdoc(mlist.matdir[m].matnum, &matval);
00269     plane=matval.plane;
00270     /* did plane change? */
00271     if(plane!=prev_plane) {seqplane++; frame=1;} else frame++;
00272     prev_plane=plane;
00273     img->planeNumber[seqplane]=plane;
00274     /* Read subheader */
00275     if(main_header.file_type==IMAGE_DATA) {
00276       ret=ecat63ReadImageheader(fp, mlist.matdir[m].strtblk, &image_header);
00277       if(dim_x!=image_header.dimension_1 || dim_y!=image_header.dimension_2) ret=1;
00278       scale=image_header.quant_scale;
00279       if(image_header.ecat_calibration_fctr>0.0
00280          && fabs(main_header.calibration_factor/image_header.ecat_calibration_fctr-1.0)>0.0001)
00281         scale*=image_header.ecat_calibration_fctr;
00282       if(img->unit==0 && image_header.quant_units>0) img->unit=image_header.quant_units;
00283       img->_dataType=image_header.data_type;
00284       img->start[frame-1]=image_header.frame_start_time/1000.;
00285       img->end[frame-1]=img->start[frame-1]+image_header.frame_duration/1000.;
00286       img->mid[frame-1]=0.5*(img->start[frame-1]+img->end[frame-1]);
00287       if(image_header.decay_corr_fctr>1.0)
00288         img->decayCorrFactor[frame-1]=image_header.decay_corr_fctr;
00289     } else if(main_header.file_type==RAW_DATA) {
00290       ret=ecat63ReadScanheader(fp, mlist.matdir[m].strtblk, &scan_header);
00291       if(dim_x!=scan_header.dimension_1 || dim_y!=scan_header.dimension_2) ret=1;
00292       scale=scan_header.scale_factor;
00293       if(scan_header.loss_correction_fctr>0.0) scale*=scan_header.loss_correction_fctr;
00294       img->_dataType=scan_header.data_type;
00295       img->start[frame-1]=scan_header.frame_start_time/1000.;
00296       img->end[frame-1]=img->start[frame-1]+scan_header.frame_duration/1000.;
00297       img->mid[frame-1]=0.5*(img->start[frame-1]+img->end[frame-1]);
00298       img->sampleDistance=10.0*scan_header.sample_distance;
00299       img->prompts[frame-1]=(float)scan_header.prompts;
00300       img->randoms[frame-1]=(float)scan_header.delayed;
00301     } else if(main_header.file_type==ATTN_DATA) {
00302       ret=ecat63ReadAttnheader(fp, mlist.matdir[m].strtblk, &attn_header);
00303       if(dim_x!=attn_header.dimension_1 || dim_y!=attn_header.dimension_2) ret=1;
00304       scale=attn_header.scale_factor;
00305       img->_dataType=attn_header.data_type;
00306     } else if(main_header.file_type==NORM_DATA) {
00307       ret=ecat63ReadNormheader(fp, mlist.matdir[m].strtblk, &norm_header);
00308       if(dim_x!=norm_header.dimension_1 || dim_y!=norm_header.dimension_2) ret=1;
00309       scale=norm_header.scale_factor;
00310       img->_dataType=norm_header.data_type;
00311     } else img->_dataType=-1;
00312     if(ret) {
00313       sprintf(ecat63errmsg, "cannot read matrix %u subheader", mlist.matdir[m].matnum);
00314       ecat63EmptyMatlist(&mlist); free(mdata); fclose(fp); return(10);
00315     }
00316     if(main_header.calibration_factor>0.0) scale*=main_header.calibration_factor;
00317     if(IMG_TEST>2) printf("scale=%e datatype=%d\n", scale, img->_dataType);
00318     /* Read the pixel data */
00319     ret=ecat63ReadMatdata(fp, mlist.matdir[m].strtblk+1,
00320          mlist.matdir[m].endblk-mlist.matdir[m].strtblk,
00321          mdata, img->_dataType);
00322     if(ret) {
00323       strcpy(ecat63errmsg, "cannot read matrix data");
00324       ecat63EmptyMatlist(&mlist); free(mdata); fclose(fp); return(11);
00325     }
00326     if(img->_dataType==BYTE_TYPE) {
00327       for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr++) {
00328         img->m[seqplane][i][j][frame-1]=scale*(float)(*mptr);
00329       }
00330     } else if(img->_dataType==VAX_I2 || img->_dataType==SUN_I2) {
00331       for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr+=2) {
00332         sptr=(short int*)mptr;
00333         img->m[seqplane][i][j][frame-1]=scale*(float)(*sptr);
00334       }
00335     } else if(img->_dataType==VAX_I4 || img->_dataType==SUN_I4) {
00336       for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr+=4) {
00337         iptr=(int*)mptr;
00338         img->m[seqplane][i][j][frame-1]=1.0; /*scale*(float)(*iptr);*/
00339       }
00340     } else if(img->_dataType==VAX_R4 || img->_dataType==IEEE_R4) {
00341       for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr+=4) {
00342         memcpy(&img->m[seqplane][i][j][frame-1], mptr, 4);
00343         img->m[seqplane][i][j][frame-1]*=scale;
00344       }
00345     }
00346   } /* next matrix */
00347 
00348   /* Set the unit */
00349   imgUnitFromEcat(img, img->unit);
00350 
00351   /* Free memory and close file */
00352   free(mdata);
00353   ecat63EmptyMatlist(&mlist);
00354   fclose(fp);
00355 
00356   /* For saving, only 2-byte VAX data types are allowed */
00357   if(img->_dataType==VAX_I4 || img->_dataType==VAX_R4) img->_dataType=VAX_I2;
00358 
00359   return(0);
00360 }
00361 /*****************************************************************************/
00362 
00363 /*****************************************************************************/
00374 int ecat63WriteAllImg(const char *fname, IMG *img) {
00375   int                 frame, plane, n, i, j, m, ret=0;
00376   float               f, fmax, fmin, g, scale;
00377   short int          *sdata, *sptr, smin, smax;
00378   FILE               *fp;
00379   ECAT63_mainheader   main_header;
00380   ECAT63_imageheader  image_header;
00381   ECAT63_scanheader   scan_header;
00382   struct tm          *scanStart;
00383 
00384 
00385   if(IMG_TEST) printf("ecat63WriteAllImg(%s, *img)\n", fname);
00386   /* Check the arguments */
00387   if(fname==NULL) {strcpy(ecat63errmsg, "invalid ECAT filename"); return(1);}
00388   if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) {
00389     strcpy(ecat63errmsg, "image data is empty"); return(2);}
00390   if(img->_dataType<1) img->_dataType=VAX_I2;
00391 
00392   /* Initiate headers */
00393   memset(&main_header, 0, sizeof(ECAT63_mainheader));
00394   memset(&image_header,0, sizeof(ECAT63_imageheader));
00395   memset(&scan_header, 0, sizeof(ECAT63_scanheader));
00396 
00397   /* Make a main header */
00398   main_header.sw_version=2;
00399   main_header.num_planes=img->dimz;  
00400   main_header.num_frames=img->dimt;  
00401   main_header.num_gates=1;
00402   main_header.num_bed_pos=1;
00403   if(img->type==IMG_TYPE_IMAGE) main_header.file_type=IMAGE_DATA;
00404   else main_header.file_type=RAW_DATA;
00405   main_header.data_type=img->_dataType;
00406   if(img->scanner>0) main_header.system_type=img->scanner;
00407   else main_header.system_type=ECAT63_SYSTEM_TYPE_DEFAULT;
00408   main_header.calibration_factor=1.0;
00409   main_header.axial_fov=img->axialFOV/10.0;
00410   main_header.transaxial_fov=img->transaxialFOV/10.0;
00411   main_header.plane_separation=img->sizez/10.0;
00412   main_header.calibration_units=imgUnitToEcat6(img);
00413   strncpy(main_header.radiopharmaceutical, img->radiopharmaceutical, 32);
00414   scanStart=localtime(&img->scanStart);
00415   if(scanStart!=NULL) {
00416     main_header.scan_start_year=scanStart->tm_year+1900;
00417     main_header.scan_start_month=scanStart->tm_mon+1;
00418     main_header.scan_start_day=scanStart->tm_mday;
00419     main_header.scan_start_hour=scanStart->tm_hour;
00420     main_header.scan_start_minute=scanStart->tm_min;
00421     main_header.scan_start_second=scanStart->tm_sec;
00422   } else {
00423     main_header.scan_start_year=1900;
00424     main_header.scan_start_month=1;
00425     main_header.scan_start_day=1;
00426     main_header.scan_start_hour=0;
00427     main_header.scan_start_minute=0;
00428     main_header.scan_start_second=0;
00429   }
00430   main_header.isotope_halflife=img->isotopeHalflife;
00431   strcpy(main_header.isotope_code, imgIsotope(img));
00432   strcpy(main_header.study_name, img->studyNr);
00433   strcpy(main_header.patient_name, img->patientName);
00434   strcpy(main_header.patient_id, img->patientID);
00435   strncpy(main_header.user_process_code, img->userProcessCode, 10);
00436   strncpy(main_header.study_description, img->studyDescription, 32);
00437   if(IMG_TEST) ecat63PrintMainheader(&main_header, stdout);
00438 
00439   /* Allocate memory for matrix data */
00440   sdata=(short int*)malloc(img->dimx*img->dimy*sizeof(short int));
00441   if(sdata==NULL) {strcpy(ecat63errmsg, "out of memory"); return(4);}
00442 
00443   /* Open output ECAT file */
00444   fp=ecat63Create(fname, &main_header);
00445   if(fp==NULL) {strcpy(ecat63errmsg, "cannot write ECAT file"); return(3);}
00446 
00447   /* Set the subheader contents, as far as possible */
00448   switch(main_header.file_type) {
00449     case RAW_DATA:
00450       scan_header.data_type=main_header.data_type;
00451       scan_header.dimension_1=img->dimx;
00452       scan_header.dimension_2=img->dimy;
00453       scan_header.frame_duration_sec=1;
00454       scan_header.scale_factor=1.0;
00455       scan_header.frame_start_time=0;
00456       scan_header.frame_duration=1000;
00457       scan_header.loss_correction_fctr=1.0;
00458       /*if(IMG_TEST) ecat63PrintScanheader(&scan_header);*/
00459       break;
00460     case IMAGE_DATA:
00461       image_header.data_type=main_header.data_type;
00462       image_header.num_dimensions=2;
00463       image_header.dimension_1=img->dimx;
00464       image_header.dimension_2=img->dimy;
00465       image_header.recon_scale=img->zoom;
00466       image_header.quant_scale=1.0;
00467       image_header.slice_width=img->sizez/10.;
00468       image_header.pixel_size=img->sizex/10.;
00469       image_header.frame_start_time=0;
00470       image_header.frame_duration=1000;
00471       image_header.plane_eff_corr_fctr=1.0;
00472       image_header.decay_corr_fctr=1.0;
00473       image_header.loss_corr_fctr=1.0;
00474       image_header.ecat_calibration_fctr=1.0;
00475       image_header.well_counter_cal_fctr=1.0;
00476       image_header.quant_units=main_header.calibration_units;
00477       /*if(IMG_TEST) ecat63PrintImageheader(&image_header);*/
00478       break;
00479   }
00480 
00481   /* Write one matrix at a time */
00482   n=0;
00483   for(plane=1; plane<=img->dimz; plane++) for(frame=1; frame<=img->dimt; frame++) {
00484     /* Scale data to short ints */
00485     /* Search min and max */
00486     fmin=fmax=f=img->m[plane-1][0][0][frame-1];
00487     for(i=0; i<img->dimy; i++) for(j=0; j<img->dimx; j++) {
00488       f=img->m[plane-1][i][j][frame-1];
00489       if(f>fmax) fmax=f; else if(f<fmin) fmin=f;
00490     }
00491     /* Calculate scaling factor */
00492     if(fabs(fmin)>fabs(fmax)) g=fabs(fmin); else g=fabs(fmax);
00493     if(g!=0) scale=32766./g; else scale=1.0;
00494     /*printf("fmin=%e fmax=%e g=%e scale=%e\n", fmin, fmax, g, scale);*/
00495     /* Scale matrix data to shorts */
00496     sptr=sdata;
00497     for(i=0; i<img->dimy; i++) for(j=0; j<img->dimx; j++) {
00498       *sptr=(short int)temp_roundf(scale*img->m[plane-1][i][j][frame-1]);
00499       sptr++;
00500     }
00501     /* Calculate and set subheader min&max and scale */
00502     smin=(short int)temp_roundf(scale*fmin); smax=(short int)temp_roundf(scale*fmax);
00503     if(main_header.file_type==RAW_DATA) {
00504       scan_header.scan_min=smin; scan_header.scan_max=smax;
00505       scan_header.scale_factor=1.0/scale;
00506       scan_header.frame_start_time=(int)temp_roundf(1000.*img->start[frame-1]);
00507       scan_header.frame_duration=
00508         (int)temp_roundf(1000.*(img->end[frame-1]-img->start[frame-1]));
00509       scan_header.sample_distance=(img->sampleDistance)/10.0;
00510       scan_header.prompts=temp_roundf(img->prompts[frame-1]);
00511       scan_header.delayed=temp_roundf(img->randoms[frame-1]);
00512     } else if(main_header.file_type==IMAGE_DATA) {
00513       image_header.image_min=smin; image_header.image_max=smax;
00514       image_header.quant_scale=1.0/scale;
00515       image_header.frame_start_time=(int)temp_roundf(1000.*img->start[frame-1]);
00516       image_header.frame_duration=
00517         (int)temp_roundf(1000.*(img->end[frame-1]-img->start[frame-1]));
00518       /* Set decay correction factor */
00519       if(img->decayCorrected)
00520         image_header.decay_corr_fctr=img->decayCorrFactor[frame-1];
00521     } 
00522     /* Write matrix data */
00523     m=mat_numcod(frame, img->planeNumber[plane-1], 1, 0, 0);
00524     /*m=mat_numcod(frame, plane, 1, 0, 0);*/
00525     sptr=sdata;
00526     if(IMG_TEST) printf("  writing matnum=%d\n", m);
00527     if(main_header.file_type==RAW_DATA) {
00528       if(IMG_TEST) ecat63PrintScanheader(&scan_header, stdout);
00529       ret=ecat63WriteScan(fp, m, &scan_header, sptr);
00530     } else if(main_header.file_type==IMAGE_DATA) {
00531       if(IMG_TEST) ecat63PrintImageheader(&image_header, stdout);
00532       ret=ecat63WriteImage(fp, m, &image_header, sptr);
00533     }
00534     if(ret) {
00535       sprintf(ecat63errmsg, "cannot write data on pl%02d fr%02d (%d).",
00536         plane, frame, ret);
00537       fclose(fp); remove(fname); free(sdata);
00538       return(9);
00539     }
00540     n++;
00541   } /* next matrix */
00542   if(IMG_TEST) printf("  %d matrices written in %s\n", n, fname);
00543 
00544   /* Close file and free memory */
00545   fclose(fp); free(sdata);
00546 
00547   return(0);
00548 }
00549 /*****************************************************************************/
00550 
00551 /*****************************************************************************/
00568 int ecat63ReadPlaneToImg(const char *fname, IMG *img) {
00569   int                 i, j, m, ret, blkNr=0, dim_x=0, dim_y=0, pxlNr=0, del_nr=0;
00570   int                 frameNr, planeNr, datatype=0, firstm=0, isFirst=0;
00571   int                 frame, plane, prev_frame, prev_plane, prev_frameNr=0;
00572   FILE               *fp;
00573   ECAT63_mainheader   main_header;
00574   MATRIXLIST          mlist;
00575   MatDir              tmpMatdir;
00576   Matval              matval;
00577   ECAT63_imageheader  image_header;
00578   ECAT63_scanheader   scan_header;
00579   ECAT63_attnheader   attn_header;
00580   ECAT63_normheader   norm_header;
00581   float               scale=1.0;
00582   short int          *sptr;
00583   char               *mdata=NULL, *mptr;
00584   int                *iptr;
00585   struct tm           scanStart;
00586 
00587 
00588   if(IMG_TEST) printf("ecat63ReadPlaneToImg(%s, *img)\n", fname);
00589   /* Check the arguments */
00590   if(fname==NULL) {strcpy(ecat63errmsg, "invalid ECAT filename"); return(2);}
00591   if(img==NULL || img->status==IMG_STATUS_UNINITIALIZED) {
00592     strcpy(ecat63errmsg, "image data not initialized"); return(2);}
00593 
00594   /* Open input CTI file for read */
00595   if((fp=fopen(fname, "rb")) == NULL) {
00596     strcpy(ecat63errmsg, "cannot open ECAT file"); return(3);}
00597 
00598   /* Read main header */
00599   if((ret=ecat63ReadMainheader(fp, &main_header))) {
00600     sprintf(ecat63errmsg, "cannot read main header (%d)", ret);
00601     fclose(fp); return(4);
00602   }
00603   if(IMG_TEST) ecat63PrintMainheader(&main_header, stdout);
00604 
00605   /* Read matrix list and nr */
00606   ecat63InitMatlist(&mlist);
00607   ret=ecat63ReadMatlist(fp, &mlist);
00608   if(ret) {
00609     sprintf(ecat63errmsg, "cannot read matrix list (%d)", ret);
00610     fclose(fp); return(5);
00611   }
00612   if(mlist.matrixNr<=0) {
00613     strcpy(ecat63errmsg, "matrix list is empty"); fclose(fp); return(5);}
00614   /* Sort matrix list */
00615   for(i=0; i<mlist.matrixNr-1; i++) for(j=i+1; j<mlist.matrixNr; j++) {
00616     if(mlist.matdir[i].matnum>mlist.matdir[j].matnum) {
00617       tmpMatdir=mlist.matdir[i];
00618       mlist.matdir[i]=mlist.matdir[j]; mlist.matdir[j]=tmpMatdir;
00619     }
00620   }
00621   /* Trim extra frames */
00622   if(main_header.num_frames>0) {
00623     del_nr=ecat63DeleteLateFrames(&mlist, main_header.num_frames);
00624     if(IMG_TEST && del_nr>0)
00625       printf("  %d entries in matrix list will not be used.\n", del_nr);
00626   }
00627 
00628   /* Decide the plane to read */
00629   if(img->status==IMG_STATUS_OCCUPIED) {
00630     /* img contains data */
00631     isFirst=0; prev_frameNr=img->dimt;
00632     /* get the current plane in there */
00633     prev_plane=img->planeNumber[img->dimz-1];
00634     /* Clear all data in img: not here but only after finding new data */
00635   } else {
00636     /* img does not contain any data */
00637     isFirst=1;
00638     /* set "current plane" to -1 */
00639     prev_plane=-1;
00640   }
00641   /* from sorted matrix list, find the first plane larger than the current one */
00642   for(m=0, plane=-1; m<mlist.matrixNr; m++) if(mlist.matdir[m].matstat==1) {
00643     mat_numdoc(mlist.matdir[m].matnum, &matval);
00644     if(matval.plane>prev_plane) {plane=matval.plane; break;}
00645   }
00646   /* If not found, return an error or 1 if this was not the first one */
00647   if(plane<0) {
00648     fclose(fp); ecat63EmptyMatlist(&mlist);
00649     if(isFirst) {
00650       sprintf(ecat63errmsg, "ECAT file contains no matrices");
00651       return(7);
00652     } else {
00653       sprintf(ecat63errmsg, "ECAT file contains no more planes");
00654       if(IMG_TEST) printf("%s\n", ecat63errmsg);
00655       return(1);
00656     }
00657   }
00658   if(IMG_TEST) printf("Next plane to read: %d\n", plane);
00659   /* Clear all data in img */
00660   imgEmpty(img);
00661 
00662   /* Calculate the number of frames and gates */
00663   prev_frame=frame=-1; frameNr=0; ret=0; planeNr=1;
00664   for(m=0; m<mlist.matrixNr; m++) if(mlist.matdir[m].matstat==1) {
00665     /* get frame and plane; work only with current plane */
00666     mat_numdoc(mlist.matdir[m].matnum, &matval);
00667     if(matval.plane<plane) continue; else if(matval.plane>plane) break;
00668     if(main_header.num_frames>=main_header.num_gates) frame=matval.frame;
00669     else frame=matval.gate;
00670     frameNr++;
00671     /* frame (gate) numbers must be continuous */
00672     if(prev_frame>0 && frame!=prev_frame+1) {ret=1; break;}
00673     prev_frame=frame;
00674     /* Calculate and check the size of one data matrix */
00675     if(frameNr==1) {
00676       firstm=m;
00677       blkNr=mlist.matdir[m].endblk-mlist.matdir[m].strtblk;
00678     } else if(blkNr!=mlist.matdir[m].endblk-mlist.matdir[m].strtblk) {
00679       ret=2; break;
00680     }
00681     prev_frame=frame;
00682   } /* next matrix */
00683   if(ret==1) {
00684     strcpy(ecat63errmsg, "missing matrix");
00685     ecat63EmptyMatlist(&mlist); fclose(fp); return(6);
00686   }
00687   if(ret==2) {
00688     strcpy(ecat63errmsg, "matrix sizes are different");
00689     ecat63EmptyMatlist(&mlist); fclose(fp); return(6);
00690   }
00691   /* Check that frameNr is equal to the dimt in IMG */
00692   if(!isFirst && frameNr!=prev_frameNr) {
00693     strcpy(ecat63errmsg, "different frame nr in different planes");
00694     ecat63EmptyMatlist(&mlist); fclose(fp); return(6);
00695   }
00696 
00697   /* Read x,y-dimensions from 1st matrix subheader on current plane */
00698   m=firstm; if(main_header.file_type==IMAGE_DATA) {
00699     ret=ecat63ReadImageheader(fp, mlist.matdir[m].strtblk, &image_header);
00700     dim_x=image_header.dimension_1; dim_y=image_header.dimension_2;
00701   } else if(main_header.file_type==RAW_DATA) {
00702     ret=ecat63ReadScanheader(fp, mlist.matdir[m].strtblk, &scan_header);
00703     dim_x=scan_header.dimension_1; dim_y=scan_header.dimension_2;
00704   } else if(main_header.file_type==ATTN_DATA) {
00705     ret=ecat63ReadAttnheader(fp, mlist.matdir[m].strtblk, &attn_header);
00706     dim_x=attn_header.dimension_1; dim_y=attn_header.dimension_2;
00707   } else if(main_header.file_type==NORM_DATA) {
00708     ret=ecat63ReadNormheader(fp, mlist.matdir[m].strtblk, &norm_header);
00709     dim_x=norm_header.dimension_1; dim_y=norm_header.dimension_2;
00710   }
00711   pxlNr=dim_x*dim_y;
00712   if(ret) {
00713     sprintf(ecat63errmsg, "cannot read matrix %u subheader", mlist.matdir[m].matnum);
00714     ecat63EmptyMatlist(&mlist); fclose(fp); return(7);
00715   }
00716 
00717   /* Allocate memory for ECAT data matrix */
00718   if(IMG_TEST) printf("allocating memory for %d blocks\n", blkNr);
00719   mdata=(char*)malloc(blkNr*MatBLKSIZE); if(mdata==NULL) {
00720     strcpy(ecat63errmsg, "out of memory");
00721     fclose(fp); ecat63EmptyMatlist(&mlist); return(8);
00722   }
00723   /* Allocate memory for whole img data */
00724   ret=imgAllocate(img, planeNr, dim_y, dim_x, frameNr);
00725   if(ret) {
00726     sprintf(ecat63errmsg, "out of memory (%d)", ret);
00727     fclose(fp); ecat63EmptyMatlist(&mlist); free(mdata); return(8);
00728   }
00729 
00730   /* Fill img info with ECAT main and sub header information */
00731   img->scanner=main_header.system_type;
00732   img->unit=main_header.calibration_units; /* this continues below */
00733   strncpy(img->radiopharmaceutical, main_header.radiopharmaceutical, 32);
00734   img->isotopeHalflife=main_header.isotope_halflife;
00735   memset(&scanStart, 0, sizeof(struct tm));
00736   scanStart.tm_year=main_header.scan_start_year-1900;
00737   scanStart.tm_mon=main_header.scan_start_month-1;
00738   scanStart.tm_mday=main_header.scan_start_day; scanStart.tm_yday=0;
00739   scanStart.tm_hour=main_header.scan_start_hour;
00740   scanStart.tm_min=main_header.scan_start_minute;
00741   scanStart.tm_sec=main_header.scan_start_second;
00742   scanStart.tm_isdst=-1;
00743   img->scanStart=mktime(&scanStart); /*printf("img->scanStart=%d\n", img->scanStart);*/
00744   if(img->scanStart==-1) img->scanStart=0;
00745   img->axialFOV=10.*main_header.axial_fov;
00746   img->transaxialFOV=10.*main_header.transaxial_fov;
00747   strncpy(img->studyNr, main_header.study_name, MAX_STUDYNR_LEN);
00748   strncpy(img->patientName, main_header.patient_name, 31);
00749   strncpy(img->patientID, main_header.patient_id, 15);
00750   img->sizez=10.*main_header.plane_separation;
00751   if(main_header.file_type==IMAGE_DATA) {
00752     img->type=IMG_TYPE_IMAGE;
00753     if(img->unit<1) img->unit=image_header.quant_units;
00754     img->zoom=image_header.recon_scale;
00755     if(image_header.decay_corr_fctr>1.0) img->decayCorrected=1;
00756     img->sizex=img->sizey=10.*image_header.pixel_size;
00757     if(img->sizez<10.*image_header.slice_width)
00758       img->sizez=10.*image_header.slice_width;
00759   } else if(main_header.file_type==RAW_DATA) {
00760     img->type=IMG_TYPE_RAW;
00761     img->decayCorrected=0;
00762   } else if(main_header.file_type==ATTN_DATA) {
00763     img->type=IMG_TYPE_RAW;
00764     img->decayCorrected=0;
00765   } else if(main_header.file_type==NORM_DATA) {
00766     img->type=IMG_TYPE_RAW;
00767     img->decayCorrected=0;
00768   }
00769   img->planeNumber[0]=plane;
00770   strncpy(img->studyDescription, main_header.study_description, 32);
00771   strncpy(img->userProcessCode, main_header.user_process_code, 10); img->userProcessCode[10]=(char)0;
00772   /* If valid study number is found in user_process_code, then take it */  
00773   if(!img->studyNr[0] && studynr_validity_check(img->userProcessCode))
00774     strcpy(img->studyNr, img->userProcessCode);
00775   /* Set _fileFormat */
00776   img->_fileFormat=IMG_E63;
00777 
00778   /* Read one ECAT matrix at a time and put them to img */
00779   for(m=firstm, frame=1; m<mlist.matrixNr; m++) if(mlist.matdir[m].matstat==1) {
00780     /* get plane */
00781     mat_numdoc(mlist.matdir[m].matnum, &matval);
00782     if(matval.plane>plane) break; /* Quit when current plane is read */
00783     /* Read subheader */
00784     if(main_header.file_type==IMAGE_DATA) {
00785       ret=ecat63ReadImageheader(fp, mlist.matdir[m].strtblk, &image_header);
00786       if(dim_x!=image_header.dimension_1 || dim_y!=image_header.dimension_2) ret=1;
00787       scale=image_header.quant_scale;
00788       if(image_header.ecat_calibration_fctr>0.0
00789          && fabs(main_header.calibration_factor/image_header.ecat_calibration_fctr-1.0)>0.0001)
00790         scale*=image_header.ecat_calibration_fctr;
00791       if(img->unit==0 && image_header.quant_units>0) img->unit=image_header.quant_units;
00792       datatype=image_header.data_type;
00793       img->start[frame-1]=image_header.frame_start_time/1000.;
00794       img->end[frame-1]=img->start[frame-1]+image_header.frame_duration/1000.;
00795       img->mid[frame-1]=0.5*(img->start[frame-1]+img->end[frame-1]);
00796       if(image_header.decay_corr_fctr>1.0)
00797         img->decayCorrFactor[frame-1]=image_header.decay_corr_fctr;
00798     } else if(main_header.file_type==RAW_DATA) {
00799       ret=ecat63ReadScanheader(fp, mlist.matdir[m].strtblk, &scan_header);
00800       if(dim_x!=scan_header.dimension_1 || dim_y!=scan_header.dimension_2) ret=1;
00801       scale=scan_header.scale_factor;
00802       if(scan_header.loss_correction_fctr>0.0) scale*=scan_header.loss_correction_fctr;
00803       datatype=scan_header.data_type;
00804       img->start[frame-1]=scan_header.frame_start_time/1000.;
00805       img->end[frame-1]=img->start[frame-1]+scan_header.frame_duration/1000.;
00806       img->mid[frame-1]=0.5*(img->start[frame-1]+img->end[frame-1]);
00807       img->sampleDistance=10.0*scan_header.sample_distance;
00808       img->prompts[frame-1]=(float)scan_header.prompts;
00809       img->randoms[frame-1]=(float)scan_header.delayed;
00810     } else if(main_header.file_type==ATTN_DATA) {
00811       ret=ecat63ReadAttnheader(fp, mlist.matdir[m].strtblk, &attn_header);
00812       if(dim_x!=attn_header.dimension_1 || dim_y!=attn_header.dimension_2) ret=1;
00813       scale=attn_header.scale_factor;
00814       datatype=attn_header.data_type;
00815     } else if(main_header.file_type==NORM_DATA) {
00816       ret=ecat63ReadNormheader(fp, mlist.matdir[m].strtblk, &norm_header);
00817       if(dim_x!=norm_header.dimension_1 || dim_y!=norm_header.dimension_2) ret=1;
00818       scale=norm_header.scale_factor;
00819       datatype=norm_header.data_type;
00820     } else datatype=-1;
00821     if(ret) {
00822       sprintf(ecat63errmsg, "cannot read matrix %u subheader", mlist.matdir[m].matnum);
00823       ecat63EmptyMatlist(&mlist); free(mdata); fclose(fp); return(7);
00824     }
00825     if(main_header.calibration_factor>0.0) scale*=main_header.calibration_factor;
00826     if(IMG_TEST>2) printf("scale=%e datatype=%d\n", scale, datatype);
00827     /* Read the pixel data */
00828     ret=ecat63ReadMatdata(fp, mlist.matdir[m].strtblk+1,
00829          mlist.matdir[m].endblk-mlist.matdir[m].strtblk,
00830          mdata, datatype);
00831     if(ret) {
00832       strcpy(ecat63errmsg, "cannot read matrix data");
00833       ecat63EmptyMatlist(&mlist); free(mdata); fclose(fp); return(9);
00834     }
00835     if(datatype==BYTE_TYPE) {
00836       for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr++) {
00837         img->m[0][i][j][frame-1]=scale*(float)(*mptr);
00838       }
00839     } else if(datatype==VAX_I2 || datatype==SUN_I2) {
00840       for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr+=2) {
00841         sptr=(short int*)mptr;
00842         img->m[0][i][j][frame-1]=scale*(float)(*sptr);
00843       }
00844     } else if(datatype==VAX_I4 || datatype==SUN_I4) {
00845       for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr+=4) {
00846         iptr=(int*)mptr;
00847         img->m[0][i][j][frame-1]=scale*(float)(*iptr);
00848       }
00849     } else if(datatype==VAX_R4 || datatype==IEEE_R4) {
00850       for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr+=4) {
00851         memcpy(&img->m[0][i][j][frame-1], mptr, 4);
00852         img->m[0][i][j][frame-1]*=scale;
00853       }
00854     }
00855     frame++;
00856   } /* next matrix (frame) */
00857   /* Set the unit */
00858   imgUnitFromEcat(img, img->unit);
00859 
00860   /* Free memory and close file */
00861   free(mdata);
00862   ecat63EmptyMatlist(&mlist);
00863   fclose(fp);
00864 
00865   /* Set _dataType if it has not been set before */
00866   if(img->_dataType<1) img->_dataType=datatype;
00867   /* For saving, only 2-byte VAX data types are allowed */
00868   if(img->_dataType==VAX_I4 || img->_dataType==VAX_R4) img->_dataType=VAX_I2;
00869 
00870   return(0);
00871 }
00872 /*****************************************************************************/
00873 
00874 /*****************************************************************************/
00886 int ecat63AddImg(const char *fname, IMG *img) {
00887   int                 n, i, j, m, ret=0, add=0;
00888   int                 frameNr, planeNr;
00889   int                 frame, plane, prev_frame, prev_plane;
00890   float               f, fmax, fmin, g, scale;
00891   short int          *sdata, *sptr, smin, smax;
00892   FILE               *fp;
00893   ECAT63_mainheader   main_header;
00894   ECAT63_imageheader  image_header;
00895   ECAT63_scanheader   scan_header;
00896   MATRIXLIST          mlist;
00897   MatDir              tmpMatdir;
00898   Matval              matval;
00899   struct tm          *scanStart;
00900 
00901 
00902   if(IMG_TEST) printf("ecat63AddImg(%s, *img)\n", fname);
00903   if(IMG_TEST>1) imgInfo(img);
00904   /* Check the arguments */
00905   if(fname==NULL) {strcpy(ecat63errmsg, "invalid ECAT filename"); return(1);}
00906   if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) {
00907     strcpy(ecat63errmsg, "image data is empty"); return(2);}
00908   if(img->_dataType<1) img->_dataType=VAX_I2;
00909 
00910   /* Initiate headers */
00911   memset(&main_header, 0, sizeof(ECAT63_mainheader));
00912   memset(&image_header,0, sizeof(ECAT63_imageheader));
00913   memset(&scan_header, 0, sizeof(ECAT63_scanheader));
00914 
00915   /* Make a main header */
00916   main_header.sw_version=2;
00917   main_header.num_planes=img->dimz;  
00918   main_header.num_frames=img->dimt;  
00919   main_header.num_gates=1;
00920   main_header.num_bed_pos=1;
00921   if(img->type==IMG_TYPE_IMAGE) main_header.file_type=IMAGE_DATA;
00922   else main_header.file_type=RAW_DATA;
00923   main_header.data_type=img->_dataType;
00924   if(img->scanner>0) main_header.system_type=img->scanner;
00925   else main_header.system_type=ECAT63_SYSTEM_TYPE_DEFAULT;
00926   main_header.calibration_factor=1.0;
00927   main_header.axial_fov=img->axialFOV/10.0;
00928   main_header.transaxial_fov=img->transaxialFOV/10.0;
00929   main_header.plane_separation=img->sizez/10.0;
00930   main_header.calibration_units=imgUnitToEcat6(img);
00931   strncpy(main_header.radiopharmaceutical, img->radiopharmaceutical, 32);
00932   scanStart=localtime(&img->scanStart);
00933   if(scanStart!=NULL) {
00934     main_header.scan_start_year=scanStart->tm_year+1900;
00935     main_header.scan_start_month=scanStart->tm_mon+1;
00936     main_header.scan_start_day=scanStart->tm_mday;
00937     main_header.scan_start_hour=scanStart->tm_hour;
00938     main_header.scan_start_minute=scanStart->tm_min;
00939     main_header.scan_start_second=scanStart->tm_sec;
00940   } else {
00941     main_header.scan_start_year=1900;
00942     main_header.scan_start_month=1;
00943     main_header.scan_start_day=1;
00944     main_header.scan_start_hour=0;
00945     main_header.scan_start_minute=0;
00946     main_header.scan_start_second=0;
00947   }
00948   main_header.isotope_halflife=img->isotopeHalflife;
00949   strcpy(main_header.isotope_code, imgIsotope(img));
00950   strcpy(main_header.study_name, img->studyNr);
00951   strcpy(main_header.patient_name, img->patientName);
00952   strcpy(main_header.patient_id, img->patientID);
00953   strncpy(main_header.study_description, img->studyDescription, 32);
00954   strncpy(main_header.user_process_code, img->userProcessCode, 10);
00955 
00956   /* Check if the ECAT file exists */
00957   if(access(fname, 0) != -1) {
00958     if(IMG_TEST) printf("Opening existing ECAT file.\n");
00959     add=1;
00960     /* Open the ECAT file */
00961     if((fp=fopen(fname, "rb+")) == NULL) {
00962       strcpy(ecat63errmsg, "cannot create ECAT file"); return(3);}
00963     /* Read main header */
00964     if((ret=ecat63ReadMainheader(fp, &main_header))) {
00965       sprintf(ecat63errmsg, "cannot read main header (%d)", ret);
00966       fclose(fp); return(3);
00967     }
00968     fflush(fp);
00969     /* Check that filetypes are matching */
00970     if((main_header.file_type==IMAGE_DATA && img->type==IMG_TYPE_IMAGE) ||
00971        (main_header.file_type==RAW_DATA && img->type==IMG_TYPE_RAW)) {
00972       /* ok */
00973     } else {
00974       sprintf(ecat63errmsg, "cannot add different filetype");
00975       fclose(fp); return(3);
00976     }
00977   } else {
00978     if(IMG_TEST) printf("ECAT file does not exist.\n");
00979     add=0;
00980     /* Create output ECAT file */
00981     fp=ecat63Create(fname, &main_header);
00982     if(fp==NULL) {strcpy(ecat63errmsg, "cannot create ECAT file"); return(3);}
00983   }
00984   if(IMG_TEST) ecat63PrintMainheader(&main_header, stdout);
00985 
00986   /* Allocate memory for matrix data */
00987   sdata=(short int*)malloc(img->dimx*img->dimy*sizeof(short int));
00988   if(sdata==NULL) {strcpy(ecat63errmsg, "out of memory"); return(4);}
00989 
00990   /* Set the subheader contents, as far as possible */
00991   switch(main_header.file_type) {
00992     case RAW_DATA:
00993       scan_header.data_type=main_header.data_type;
00994       scan_header.dimension_1=img->dimx;
00995       scan_header.dimension_2=img->dimy;
00996       scan_header.frame_duration_sec=1;
00997       scan_header.scale_factor=1.0;
00998       scan_header.frame_start_time=0;
00999       scan_header.frame_duration=1000;
01000       scan_header.loss_correction_fctr=1.0;
01001       scan_header.sample_distance=(img->sampleDistance)/10.0;
01002       /*if(IMG_TEST) ecat63PrintScanheader(&scan_header);*/
01003       break;
01004     case IMAGE_DATA:
01005       image_header.data_type=main_header.data_type;
01006       image_header.num_dimensions=2;
01007       image_header.dimension_1=img->dimx;
01008       image_header.dimension_2=img->dimy;
01009       image_header.recon_scale=img->zoom;
01010       image_header.quant_scale=1.0;
01011       image_header.slice_width=img->sizez/10.;
01012       image_header.pixel_size=img->sizex/10.;
01013       image_header.frame_start_time=0;
01014       image_header.frame_duration=1000;
01015       image_header.plane_eff_corr_fctr=1.0;
01016       image_header.decay_corr_fctr=0.0;
01017       image_header.loss_corr_fctr=1.0;
01018       image_header.ecat_calibration_fctr=1.0;
01019       image_header.well_counter_cal_fctr=1.0;
01020       image_header.quant_units=main_header.calibration_units;
01021       /*if(IMG_TEST) ecat63PrintImageheader(&image_header);*/
01022       break;
01023   }
01024 
01025   /* Write one matrix at a time */
01026   n=0;
01027   for(plane=1; plane<=img->dimz; plane++) for(frame=1; frame<=img->dimt; frame++) {
01028     /* Scale data to short ints */
01029     /* Search min and max */
01030     fmin=fmax=f=img->m[plane-1][0][0][frame-1];
01031     for(i=0; i<img->dimy; i++) for(j=0; j<img->dimx; j++) {
01032       f=img->m[plane-1][i][j][frame-1];
01033       if(f>fmax) fmax=f; else if(f<fmin) fmin=f;
01034     }
01035     /* Calculate scaling factor */
01036     if(fabs(fmin)>fabs(fmax)) g=fabs(fmin); else g=fabs(fmax);
01037     if(g!=0) scale=32766./g; else scale=1.0;
01038     /*printf("fmin=%e fmax=%e g=%e scale=%e\n", fmin, fmax, g, scale);*/
01039     /* Scale matrix data to shorts */
01040     sptr=sdata;
01041     for(i=0; i<img->dimy; i++) for(j=0; j<img->dimx; j++) {
01042       *sptr=(short int)temp_roundf(scale*img->m[plane-1][i][j][frame-1]);
01043       sptr++;
01044     }
01045     /* Calculate and set subheader min&max and scale */
01046     smin=(short int)temp_roundf(scale*fmin); smax=(short int)temp_roundf(scale*fmax);
01047     if(main_header.file_type==RAW_DATA) {
01048       scan_header.scan_min=smin; scan_header.scan_max=smax;
01049       scan_header.scale_factor=1.0/scale;
01050       scan_header.frame_start_time=(int)temp_roundf(1000.*img->start[frame-1]);
01051       scan_header.frame_duration=
01052         (int)temp_roundf(1000.*(img->end[frame-1]-img->start[frame-1]));
01053       scan_header.prompts=temp_roundf(img->prompts[frame-1]);
01054       scan_header.delayed=temp_roundf(img->randoms[frame-1]);
01055     } else if(main_header.file_type==IMAGE_DATA) {
01056       image_header.image_min=smin; image_header.image_max=smax;
01057       image_header.quant_scale=1.0/scale;
01058       image_header.frame_start_time=(int)temp_roundf(1000.*img->start[frame-1]);
01059       image_header.frame_duration=
01060         (int)temp_roundf(1000.*(img->end[frame-1]-img->start[frame-1]));
01061       /* Set decay correction factor */
01062       if(img->decayCorrected)
01063         image_header.decay_corr_fctr=img->decayCorrFactor[frame-1];
01064     } 
01065     /* Write matrix data */
01066     m=mat_numcod(frame, img->planeNumber[plane-1], 1, 0, 0);
01067     sptr=sdata;
01068     if(IMG_TEST) printf("  writing matnum=%d\n", m);
01069     if(main_header.file_type==RAW_DATA) {
01070       /*if(IMG_TEST) ecat63PrintScanheader(&scan_header);*/
01071       ret=ecat63WriteScan(fp, m, &scan_header, sptr);
01072     } else if(main_header.file_type==IMAGE_DATA) {
01073       /*if(IMG_TEST) ecat63PrintImageheader(&image_header);*/
01074       ret=ecat63WriteImage(fp, m, &image_header, sptr);
01075     }
01076     if(ret) {
01077       sprintf(ecat63errmsg, "cannot write data on pl%02d fr%02d (%d).",
01078         plane, frame, ret);
01079       fclose(fp); remove(fname); free(sdata);
01080       return(9);
01081     }
01082     n++;
01083   } /* next matrix */
01084   if(IMG_TEST) printf("  %d matrices written in %s\n", n, fname);
01085 
01086   /* Free matrix memory */
01087   free(sdata);
01088 
01089   /* If matrices were added, set main header */
01090   if(add==1) {
01091     fflush(fp);
01092     /* Read matrix list */
01093     ecat63InitMatlist(&mlist);
01094     ret=ecat63ReadMatlist(fp, &mlist);
01095     if(ret) {
01096       sprintf(ecat63errmsg, "cannot read matrix list (%d)", ret);
01097       fclose(fp); return(21);
01098     }
01099     if(mlist.matrixNr<=0) {
01100       strcpy(ecat63errmsg, "matrix list is empty"); fclose(fp); return(21);}
01101     /* Sort matrix list */
01102     for(i=0; i<mlist.matrixNr-1; i++) for(j=i+1; j<mlist.matrixNr; j++) {
01103       if(mlist.matdir[i].matnum>mlist.matdir[j].matnum) {
01104         tmpMatdir=mlist.matdir[i];
01105         mlist.matdir[i]=mlist.matdir[j]; mlist.matdir[j]=tmpMatdir;
01106       }
01107     }
01108     /* Calculate the number of planes and frames */
01109     prev_plane=plane=-1; prev_frame=frame=-1; frameNr=planeNr=0;
01110     for(m=0; m<mlist.matrixNr; m++) {
01111       mat_numdoc(mlist.matdir[m].matnum, &matval);
01112       plane=matval.plane; frame=matval.frame;
01113       if(plane!=prev_plane) {frameNr=1; planeNr++;} else {frameNr++;}
01114       prev_plane=plane; prev_frame=frame;
01115     } /* next matrix */
01116     ecat63EmptyMatlist(&mlist);
01117     main_header.num_planes=planeNr; main_header.num_frames=frameNr;
01118     /* Write main header */
01119     ret=ecat63WriteMainheader(fp, &main_header);
01120     if(ret) {
01121       sprintf(ecat63errmsg, "cannot write mainheader (%d)", ret);
01122       fclose(fp); return(22);
01123     }
01124   }
01125 
01126   /* Close file and free memory */
01127   fclose(fp);
01128 
01129   return(0);
01130 }
01131 /*****************************************************************************/
01132 
01133 /*****************************************************************************/
01140 int imgEcat63Supported(ECAT63_mainheader *h) {
01141   if(h==NULL) return(0);
01142   if(h->file_type==IMAGE_DATA) return(1);
01143   if(h->file_type==RAW_DATA)   return(1);
01144   if(h->file_type==ATTN_DATA)  return(1);
01145   if(h->file_type==NORM_DATA)  return(1);
01146   return(0);
01147 }
01148 /*****************************************************************************/
01149 
01150 /*****************************************************************************/
01157 void imgGetEcat63MHeader(IMG *img, ECAT63_mainheader *h) {
01158   struct tm scanStart;
01159 
01160   img->_dataType=h->data_type; /* again in subheaders*/
01161   /* For saving IMG data, only 2-byte VAX data types are allowed, so change it now */
01162   if(img->_dataType==VAX_I4 || img->_dataType==VAX_R4) img->_dataType=VAX_I2;
01163   img->scanner=h->system_type;
01164   strncpy(img->radiopharmaceutical, h->radiopharmaceutical, 32);
01165   img->isotopeHalflife=h->isotope_halflife;
01166   {
01167   memset(&scanStart, 0, sizeof(struct tm));
01168   scanStart.tm_year=h->scan_start_year-1900;
01169   scanStart.tm_mon=h->scan_start_month-1;
01170   scanStart.tm_mday=h->scan_start_day; scanStart.tm_yday=0;
01171   scanStart.tm_hour=h->scan_start_hour;
01172   scanStart.tm_min=h->scan_start_minute;
01173   scanStart.tm_sec=h->scan_start_second;
01174   scanStart.tm_isdst=-1;
01175   img->scanStart=mktime(&scanStart); if(img->scanStart==-1) img->scanStart=0;
01176   }
01177   img->axialFOV=10.0*h->axial_fov;
01178   img->transaxialFOV=10.0*h->transaxial_fov;
01179   strncpy(img->studyNr, h->study_name, MAX_STUDYNR_LEN);
01180   strncpy(img->patientName, h->patient_name, 31);
01181   strncpy(img->patientID, h->patient_id, 15);
01182   img->sizez=10.0*h->plane_separation;
01183   switch(h->file_type) {
01184     case IMAGE_DATA: img->type=IMG_TYPE_IMAGE; break;
01185     case RAW_DATA:
01186     case ATTN_DATA:
01187     case NORM_DATA: 
01188     default: img->type=IMG_TYPE_RAW;
01189   }
01190   strncpy(img->studyDescription, h->study_description, 32);
01191   strncpy(img->userProcessCode, h->user_process_code, 10);
01192   img->userProcessCode[10]=(char)0;
01193   /* If valid study number is found in user_process_code, then take it */  
01194   if(!img->studyNr[0] && studynr_validity_check(img->userProcessCode))
01195     strcpy(img->studyNr, img->userProcessCode);
01196   /* Set calibration unit (this may have to be read from subheader later) */
01197   imgUnitFromEcat(img, h->calibration_units);
01198 }
01199 /*****************************************************************************/
01200 
01201 /*****************************************************************************/
01208 void imgSetEcat63MHeader(IMG *img, ECAT63_mainheader *h) {
01209   struct tm *scanStart;
01210 
01211   h->sw_version=2;
01212   h->num_planes=img->dimz;  
01213   h->num_frames=img->dimt;  
01214   h->num_gates=1;
01215   h->num_bed_pos=1;
01216   if(img->type==IMG_TYPE_IMAGE) h->file_type=IMAGE_DATA;
01217   else h->file_type=RAW_DATA;
01218   h->data_type=VAX_I2;
01219   if(img->scanner>0) h->system_type=img->scanner;
01220   else h->system_type=ECAT63_SYSTEM_TYPE_DEFAULT;
01221   h->calibration_factor=1.0;
01222   h->axial_fov=img->axialFOV/10.0;
01223   h->transaxial_fov=img->transaxialFOV/10.0;
01224   h->plane_separation=img->sizez/10.0;
01225   h->calibration_units=imgUnitToEcat6(img);
01226   strncpy(h->radiopharmaceutical, img->radiopharmaceutical, 32);
01227   scanStart=localtime(&img->scanStart);
01228   if(scanStart!=NULL) {
01229     h->scan_start_year=scanStart->tm_year+1900;
01230     h->scan_start_month=scanStart->tm_mon+1;
01231     h->scan_start_day=scanStart->tm_mday;
01232     h->scan_start_hour=scanStart->tm_hour;
01233     h->scan_start_minute=scanStart->tm_min;
01234     h->scan_start_second=scanStart->tm_sec;
01235   } else {
01236     h->scan_start_year=1900;
01237     h->scan_start_month=1;
01238     h->scan_start_day=1;
01239     h->scan_start_hour=0;
01240     h->scan_start_minute=0;
01241     h->scan_start_second=0;
01242   }
01243   h->isotope_halflife=img->isotopeHalflife;
01244   strcpy(h->isotope_code, imgIsotope(img));
01245   strcpy(h->study_name, img->studyNr);
01246   strcpy(h->patient_name, img->patientName);
01247   strcpy(h->patient_id, img->patientID);
01248   strncpy(h->user_process_code, img->userProcessCode, 10);
01249   strncpy(h->study_description, img->studyDescription, 32);
01250 }
01251 /*****************************************************************************/
01252 
01253 /*****************************************************************************/
01260 int imgGetEcat63Fileformat(ECAT63_mainheader *h) {
01261   int fileFormat=IMG_UNKNOWN;
01262   switch(h->file_type) {
01263     case IMAGE_DATA:
01264     case RAW_DATA:
01265     case ATTN_DATA:
01266     case NORM_DATA:
01267       fileFormat=IMG_E63; break;
01268     default:
01269       fileFormat=IMG_UNKNOWN; break;
01270   }
01271   return fileFormat;
01272 }
01273 /*****************************************************************************/
01274 
01275 /*****************************************************************************/
01289 int imgReadEcat63Header(const char *fname, IMG *img) {
01290   int                 m, blkNr=0, ret=0;
01291   int                 frameNr, planeNr, del_nr=0;
01292   FILE               *fp;
01293   ECAT63_mainheader   main_header;
01294   MATRIXLIST          mlist;
01295   ECAT63_imageheader  image_header;
01296   ECAT63_scanheader   scan_header;
01297   ECAT63_attnheader   attn_header;
01298   ECAT63_normheader   norm_header;
01299 
01300 
01301   if(IMG_TEST) printf("\nimgReadEcat63Header(%s, *img)\n", fname);
01302   
01303   /* Check the arguments */
01304   if(img==NULL) return STATUS_FAULT;
01305   if(img->status!=IMG_STATUS_INITIALIZED) return STATUS_FAULT;
01306   imgSetStatus(img, STATUS_FAULT);
01307   if(fname==NULL) return STATUS_FAULT;
01308     
01309   /* Open the file */
01310   if((fp=fopen(fname, "rb")) == NULL) return STATUS_NOFILE;
01311   
01312   /* Read main header */
01313   if((ret=ecat63ReadMainheader(fp, &main_header))) {
01314     fclose(fp); return STATUS_NOMAINHEADER;}
01315   /* Check if file_type is supported */
01316   if(imgEcat63Supported(&main_header)==0) {fclose(fp); return STATUS_UNSUPPORTED;}
01317   /* Copy main header information into IMG; sets also img.type */
01318   imgGetEcat63MHeader(img, &main_header);
01319   if(IMG_TEST>7) printf("img.type := %d\n", img->type);
01320   img->_fileFormat=imgGetEcat63Fileformat(&main_header);
01321   if(IMG_TEST>7) printf("img._fileFormat := %d\n", img->_fileFormat);
01322   if(img->_fileFormat==IMG_UNKNOWN) {fclose(fp); return STATUS_UNSUPPORTED;}
01323 
01324   /* Read matrix list and nr and sort it */
01325   ecat63InitMatlist(&mlist);
01326   ret=ecat63ReadMatlist(fp, &mlist);
01327   if(ret) {fclose(fp); return STATUS_NOMATLIST;}
01328   if(mlist.matrixNr<=0) {
01329     ecat63EmptyMatlist(&mlist); fclose(fp); return STATUS_INVALIDMATLIST;}
01330   /* Make sure that plane and frame numbers are continuous */
01331   ecat63GatherMatlist(&mlist, 1, 1, 1, 1);
01332   ecat63SortMatlistByPlane(&mlist); /* otherwise frameNr cannot be computed as below */
01333   /* Trim extra frames */
01334   if(main_header.num_frames>0)
01335     del_nr=ecat63DeleteLateFrames(&mlist, main_header.num_frames);
01336   /* Get plane and frame numbers and check that volume is full */
01337   ret=ecat63GetPlaneAndFrameNr(&mlist, &main_header, &planeNr, &frameNr);
01338   if(ret) {ecat63EmptyMatlist(&mlist); fclose(fp); return ret;}
01339   img->dimz=planeNr;
01340   img->dimt=frameNr;
01341   /* Calculate and check the size of one data matrix */
01342   ret=ecat63GetMatrixBlockSize(&mlist, &blkNr);
01343   if(ret) {ecat63EmptyMatlist(&mlist); fclose(fp); return ret;}
01344 
01345   /* Read one subheader */
01346   if(IMG_TEST>5) printf("main_header.file_type := %d\n", main_header.file_type);
01347   m=0;
01348   switch(main_header.file_type) {
01349     case IMAGE_DATA:
01350       ret=ecat63ReadImageheader(fp, mlist.matdir[m].strtblk, &image_header);
01351       break;
01352     case RAW_DATA:
01353       ret=ecat63ReadScanheader(fp, mlist.matdir[m].strtblk, &scan_header);
01354       break;
01355     case ATTN_DATA:
01356       ret=ecat63ReadAttnheader(fp, mlist.matdir[m].strtblk, &attn_header);
01357       break;
01358     case NORM_DATA:
01359       ret=ecat63ReadNormheader(fp, mlist.matdir[m].strtblk, &norm_header);
01360       break;
01361     default: ret=-1;
01362   }
01363   /* Free locally allocated memory and close the file */
01364   ecat63EmptyMatlist(&mlist); fclose(fp);
01365   /* Check whether subheader was read */
01366   if(ret) return STATUS_NOSUBHEADER;
01367 
01368   /* Get the following information in the subheader:
01369      dimensions x, y and z; datatype;
01370      image decay correction on/off, zoom, and pixel size;
01371      sinogram sample distance.
01372    */
01373   switch(main_header.file_type) {
01374     case IMAGE_DATA:
01375       img->dimx=image_header.dimension_1; img->dimy=image_header.dimension_2;
01376       if(img->unit<1) imgUnitFromEcat(img, image_header.quant_units);
01377       img->_dataType=image_header.data_type;
01378       img->zoom=image_header.recon_scale;
01379       if(image_header.decay_corr_fctr>1.0) img->decayCorrected=1;
01380       img->sizex=img->sizey=10.*image_header.pixel_size;
01381       if(img->sizez<10.*image_header.slice_width)
01382         img->sizez=10.*image_header.slice_width;
01383       break;
01384     case RAW_DATA:
01385       img->dimx=scan_header.dimension_1; img->dimy=scan_header.dimension_2;
01386       img->type=IMG_TYPE_RAW;
01387       img->_dataType=scan_header.data_type;
01388       img->decayCorrected=0;
01389       img->sampleDistance=10.0*scan_header.sample_distance;
01390       break;
01391     case ATTN_DATA:
01392       img->dimx=attn_header.dimension_1; img->dimy=attn_header.dimension_2;
01393       img->type=IMG_TYPE_RAW;
01394       img->decayCorrected=0;
01395       img->_dataType=attn_header.data_type;
01396       break;
01397     case NORM_DATA:
01398       img->dimx=norm_header.dimension_1; img->dimy=norm_header.dimension_2;
01399       img->type=IMG_TYPE_RAW;
01400       img->decayCorrected=0;
01401       img->_dataType=norm_header.data_type;
01402       break;
01403   }
01404 
01405   /* For saving IMG data, only 2-byte VAX data types are allowed, so change it now */
01406   if(img->_dataType==VAX_I4 || img->_dataType==VAX_R4) img->_dataType=VAX_I2;
01407 
01408   imgSetStatus(img, STATUS_OK);
01409   return STATUS_OK;
01410 }
01411 /*****************************************************************************/
01412 
01413 /*****************************************************************************/
01422 int imgReadEcat63FirstFrame(const char *fname, IMG *img) {
01423   int ret=0;
01424 
01425   if(IMG_TEST) printf("\nimgReadEcat63FirstFrame(%s, *img)\n", fname);
01426   /* Check the input */
01427   if(img==NULL) return STATUS_FAULT;
01428   if(img->status!=IMG_STATUS_INITIALIZED) return STATUS_FAULT;
01429   imgSetStatus(img, STATUS_FAULT);
01430   if(fname==NULL) return STATUS_FAULT;
01431 
01432   /* Read header information from file */
01433   ret=imgReadEcat63Header(fname, img); if(ret) return(ret);
01434   if(IMG_TEST>3) imgInfo(img);
01435 
01436   /* Allocate memory for one frame */
01437   img->dimt=1;
01438   ret=imgAllocate(img, img->dimz, img->dimy, img->dimx, img->dimt);
01439   if(ret) return STATUS_NOMEMORY;
01440 
01441   /* Read the first frame */
01442   ret=imgReadEcat63Frame(fname, 1, img, 0); if(ret) return(ret); 
01443 
01444   /* All went well */
01445   imgSetStatus(img, STATUS_OK);
01446   return STATUS_OK;
01447 }
01448 /*****************************************************************************/
01449 
01450 /*****************************************************************************/
01464 int imgReadEcat63Frame(const char *fname, int frame_to_read, IMG *img, int frame_index) {
01465   int                 m, ret=0, blkNr=0, frame, plane, seqplane=-1, xi, yi;
01466   int                 local_data_type=0;
01467   FILE               *fp;
01468   ECAT63_mainheader   main_header;
01469   MATRIXLIST          mlist;
01470   Matval              matval;
01471   ECAT63_imageheader  image_header;
01472   ECAT63_scanheader   scan_header;
01473   ECAT63_attnheader   attn_header;
01474   ECAT63_normheader   norm_header;
01475   float               scale=1.0;
01476   short int          *sptr;
01477   char               *mdata=NULL, *mptr;
01478   int                *iptr;
01479 
01480 
01481   if(IMG_TEST) printf("\nimgReadEcat63Frame(%s, %d, *img, %d)\n",
01482     fname, frame_to_read, frame_index);
01483     
01484   /* Check the input */
01485   if(img==NULL) return STATUS_FAULT;
01486   if(img->status!=IMG_STATUS_OCCUPIED) return STATUS_FAULT;
01487   imgSetStatus(img, STATUS_FAULT);
01488   if(fname==NULL) return STATUS_FAULT;
01489   if(frame_index<0 || frame_index>img->dimt-1) return STATUS_FAULT;
01490   if(frame_to_read<1) return STATUS_FAULT;
01491 
01492   /* Open the file */
01493   if((fp=fopen(fname, "rb")) == NULL) return STATUS_NOFILE;
01494   
01495   /* Read main header */
01496   if((ret=ecat63ReadMainheader(fp, &main_header))) {
01497     fclose(fp); return STATUS_NOMAINHEADER;}
01498 
01499   /* Read matrix list and nr and sort it */
01500   ecat63InitMatlist(&mlist);
01501   ret=ecat63ReadMatlist(fp, &mlist);
01502   if(ret) {fclose(fp); return STATUS_NOMATLIST;}
01503   if(mlist.matrixNr<=0) {
01504     fclose(fp); ecat63EmptyMatlist(&mlist); return STATUS_INVALIDMATLIST;}
01505   /* Make sure that plane and frame numbers are continuous */
01506   ecat63GatherMatlist(&mlist, 1, 1, 1, 1);
01507   ecat63SortMatlistByFrame(&mlist);
01508   
01509   /* Calculate and check the size of one data matrix */
01510   ret=ecat63GetMatrixBlockSize(&mlist, &blkNr);
01511   if(ret) {ecat63EmptyMatlist(&mlist); fclose(fp); return ret;}
01512   /* And allocate memory for the raw data blocks */
01513   if(IMG_TEST>6) printf("allocating memory for %d blocks\n", blkNr);
01514   mdata=(char*)malloc(blkNr*MatBLKSIZE);
01515   if(mdata==NULL) {fclose(fp); ecat63EmptyMatlist(&mlist); return STATUS_NOMEMORY;}
01516 
01517   /* Read all matrices that belong to the required frame */
01518   ret=0; seqplane=-1;
01519   for(m=0; m<mlist.matrixNr; m++) if(mlist.matdir[m].matstat==1) {
01520     /* get frame and plane */
01521     mat_numdoc(mlist.matdir[m].matnum, &matval);
01522     plane=matval.plane;
01523     if(main_header.num_frames>=main_header.num_gates) frame=matval.frame;
01524     else frame=matval.gate; /* printf("frame=%d plane=%d\n", frame, plane); */
01525     if(frame!=frame_to_read) continue; /* do not process other frames */
01526     seqplane++; 
01527     if(img->planeNumber[seqplane]<1) {
01528       img->planeNumber[seqplane]=plane;
01529     } else if(img->planeNumber[seqplane]!=plane) {
01530       fclose(fp); ecat63EmptyMatlist(&mlist); free(mdata);
01531       return STATUS_MISSINGMATRIX;    
01532     }
01533     
01534     /* Read the subheader */
01535     if(IMG_TEST>4) printf("reading subheader for matrix %d,%d\n", frame, plane);
01536     if(main_header.file_type==IMAGE_DATA) {
01537       ret=ecat63ReadImageheader(fp, mlist.matdir[m].strtblk, &image_header);
01538       scale=image_header.quant_scale;
01539       if(image_header.ecat_calibration_fctr>0.0
01540          && fabs(main_header.calibration_factor/image_header.ecat_calibration_fctr-1.0)>0.0001)
01541         scale*=image_header.ecat_calibration_fctr;
01542       local_data_type=image_header.data_type;
01543       img->start[frame_index]=image_header.frame_start_time/1000.;
01544       img->end[frame_index]=img->start[frame_index]+image_header.frame_duration/1000.;
01545       img->mid[frame_index]=0.5*(img->start[frame_index]+img->end[frame_index]);
01546       if(image_header.decay_corr_fctr>1.0)
01547         img->decayCorrFactor[frame_index]=image_header.decay_corr_fctr;
01548     } else if(main_header.file_type==RAW_DATA) {
01549       ret=ecat63ReadScanheader(fp, mlist.matdir[m].strtblk, &scan_header);
01550       scale=scan_header.scale_factor;
01551       if(scan_header.loss_correction_fctr>0.0) scale*=scan_header.loss_correction_fctr;
01552       local_data_type=scan_header.data_type;
01553       img->start[frame_index]=scan_header.frame_start_time/1000.;
01554       img->end[frame_index]=img->start[frame_index]+scan_header.frame_duration/1000.;
01555       img->mid[frame_index]=0.5*(img->start[frame_index]+img->end[frame_index]);
01556       img->sampleDistance=10.0*scan_header.sample_distance;
01557       img->prompts[frame_index]=(float)scan_header.prompts;
01558       img->randoms[frame_index]=(float)scan_header.delayed;
01559     } else if(main_header.file_type==ATTN_DATA) {
01560       ret=ecat63ReadAttnheader(fp, mlist.matdir[m].strtblk, &attn_header);
01561       scale=attn_header.scale_factor;
01562       local_data_type=attn_header.data_type;
01563     } else if(main_header.file_type==NORM_DATA) {
01564       ret=ecat63ReadNormheader(fp, mlist.matdir[m].strtblk, &norm_header);
01565       scale=norm_header.scale_factor;
01566       local_data_type=norm_header.data_type;
01567     } else
01568       img->_dataType=-1;
01569     if(ret) {
01570       ecat63EmptyMatlist(&mlist); free(mdata); fclose(fp);
01571       return STATUS_NOSUBHEADER;
01572     }
01573     if(main_header.calibration_factor>0.0) scale*=main_header.calibration_factor;
01574     if(IMG_TEST>6) printf("scale=%e datatype=%d\n", scale, img->_dataType);
01575     /* Read the pixel data */
01576     if(IMG_TEST>4) printf("reading matrix data\n");
01577     ret=ecat63ReadMatdata(fp, mlist.matdir[m].strtblk+1,
01578          mlist.matdir[m].endblk-mlist.matdir[m].strtblk,
01579          mdata, local_data_type);
01580     if(ret) {
01581       ecat63EmptyMatlist(&mlist); free(mdata); fclose(fp);
01582       return STATUS_MISSINGMATRIX;
01583     }
01584     if(IMG_TEST>4) printf("converting matrix data to floats\n");
01585     if(local_data_type==BYTE_TYPE) {
01586       for(yi=0, mptr=mdata; yi<img->dimy; yi++)
01587         for(xi=0; xi<img->dimx; xi++, mptr++) {
01588           img->m[seqplane][yi][xi][frame_index]=scale*(float)(*mptr);
01589         }
01590     } else if(local_data_type==VAX_I2 || local_data_type==SUN_I2) {
01591       for(yi=0, mptr=mdata; yi<img->dimy; yi++)
01592         for(xi=0; xi<img->dimx; xi++, mptr+=2) {
01593           sptr=(short int*)mptr;
01594           img->m[seqplane][yi][xi][frame_index]=scale*(float)(*sptr);
01595         }
01596     } else if(local_data_type==VAX_I4 || local_data_type==SUN_I4) {
01597       for(yi=0, mptr=mdata; yi<img->dimy; yi++)
01598         for(xi=0; xi<img->dimx; xi++, mptr+=4) {
01599           iptr=(int*)mptr;
01600           img->m[seqplane][yi][xi][frame_index]=1.0; /* scale*(float)(*iptr); */
01601         }
01602     } else if(local_data_type==VAX_R4 || local_data_type==IEEE_R4) {
01603       for(yi=0, mptr=mdata; yi<img->dimy; yi++)
01604         for(xi=0; xi<img->dimx; xi++, mptr+=4) {
01605           memcpy(&img->m[seqplane][yi][xi][frame_index], mptr, 4);
01606           img->m[seqplane][yi][xi][frame_index]*=scale;
01607         }
01608     }
01609   } /* next matrix */
01610   if(IMG_TEST>3) printf("end of matrices.\n");
01611   
01612   free(mdata);
01613   ecat63EmptyMatlist(&mlist); 
01614   fclose(fp);
01615 
01616   /* seqplane is <0 only if this frame did not exist at all; return -1 */
01617   if(IMG_TEST>4) printf("last_seqplane := %d.\n", seqplane);
01618   if(seqplane<0) {imgSetStatus(img, STATUS_NOMATRIX); return STATUS_NOMATRIX;}
01619 
01620   /* check that correct number of planes was read */
01621   if(seqplane+1 != img->dimz) {
01622     imgSetStatus(img, STATUS_MISSINGMATRIX);
01623     return STATUS_MISSINGMATRIX;
01624   }
01625 
01626   /* For saving IMG data, only 2-byte VAX data types are allowed, so change it now */
01627   if(img->_dataType==VAX_I4 || img->_dataType==VAX_R4) img->_dataType=VAX_I2;
01628 
01629   /* All went well */
01630   imgSetStatus(img, STATUS_OK);
01631   return STATUS_OK;
01632 }
01633 /*****************************************************************************/
01634 
01635 /*****************************************************************************/
01656 int imgWriteEcat63Frame(const char *fname, int frame_to_write, IMG *img, int frame_index) {
01657   IMG test_img;
01658   int ret=0, pxlNr, zi, xi, yi, matrixId;
01659   ECAT63_mainheader   main_header;
01660   ECAT63_imageheader  image_header;
01661   ECAT63_scanheader   scan_header;
01662   void *sub_header=NULL;
01663   FILE *fp;
01664   float *fdata=NULL, *fptr;
01665 
01666 
01667   if(IMG_TEST) printf("\nimgWriteEcat63Frame(%s, %d, *img, %d)\n",
01668     fname, frame_to_write, frame_index);
01669   if(IMG_TEST>1) ECAT63_TEST=IMG_TEST-1;
01670     
01671   /*
01672    *  Check the input 
01673    */
01674   if(fname==NULL) return STATUS_FAULT;
01675   if(img==NULL) return STATUS_FAULT;
01676   if(img->status!=IMG_STATUS_OCCUPIED) return STATUS_FAULT;
01677   if(frame_to_write<0) return STATUS_FAULT;
01678   if(frame_index<0 || frame_index>=img->dimt) return STATUS_FAULT;
01679   if(img->_fileFormat!=IMG_E63) return STATUS_FAULT;
01680 
01681   /* Initiate headers */
01682   memset(&main_header, 0, sizeof(ECAT63_mainheader));
01683   memset(&image_header,0, sizeof(ECAT63_imageheader));
01684   memset(&scan_header, 0, sizeof(ECAT63_scanheader));
01685   imgInit(&test_img);
01686 
01687 
01688   /*
01689    *  If file does not exist, then create it with new mainheader,
01690    *  and if it does exist, then read and check header information.
01691    *  Create or edit mainheader to contain correct frame nr.   
01692    *  In any case, leave file pointer open for write.   
01693    */
01694   if(access(fname, 0) == -1) { /* file does not exist*/
01695 
01696     /* Set main header */
01697     imgSetEcat63MHeader(img, &main_header);
01698     if(IMG_TEST>2) {
01699       ecat63PrintMainheader(&main_header, stdout);
01700     }
01701     if(frame_to_write==0) frame_to_write=1;
01702     main_header.num_frames=frame_to_write;
01703 
01704     /* Open file, write main header and initiate matrix list */
01705     fp=ecat63Create(fname, &main_header); if(fp==NULL) return STATUS_NOWRITEPERM;
01706 
01707   } else { /* file exists*/
01708   
01709     /* Read header information for checking */
01710     ret=imgReadEcat63Header(fname, &test_img); if(ret!=0) return ret;
01711     /* Check that file format is the same */
01712     if(img->_fileFormat!=test_img._fileFormat || img->type!=test_img.type)
01713       return STATUS_WRONGFILETYPE;
01714     /* Check that matrix sizes are the same */
01715     if(img->dimz!=test_img.dimz || img->dimx!=test_img.dimx ||
01716        img->dimy!=test_img.dimy)
01717       return STATUS_VARMATSIZE;
01718     imgEmpty(&test_img);
01719 
01720     /* Open the file for read and write */
01721     if((fp=fopen(fname, "r+b"))==NULL) return STATUS_NOWRITEPERM;
01722 
01723     /* Read the mainheader, set new frame number, and write it back */
01724     if((ret=ecat63ReadMainheader(fp, &main_header))!=0) return STATUS_NOMAINHEADER;
01725     if(frame_to_write==0) frame_to_write=main_header.num_frames+1;
01726     if(main_header.num_frames<frame_to_write)
01727       main_header.num_frames=frame_to_write;
01728     if((ret=ecat63WriteMainheader(fp, &main_header))!=0) return STATUS_NOWRITEPERM;
01729     
01730   }
01731   if(IMG_TEST>2) {
01732     printf("frame_to_write := %d\n", frame_to_write);
01733   }
01734 
01735   /* Allocate memory for matrix float data */
01736   pxlNr=img->dimx*img->dimy*img->dimz;
01737   fdata=(float*)malloc(pxlNr*sizeof(float));
01738   if(fdata==NULL) {fclose(fp); return STATUS_NOMEMORY;}
01739   
01740   /* Set matrix subheader */
01741   if(img->type==IMG_TYPE_RAW) sub_header=(void*)&scan_header;
01742   else if(img->type==IMG_TYPE_IMAGE) sub_header=&image_header;
01743   else {fclose(fp); return STATUS_FAULT;}
01744   imgSetEcat63SHeader(img, sub_header);
01745 
01746   /* Copy matrix pixel values to fdata */
01747   fptr=fdata;
01748   for(zi=0; zi<img->dimz; zi++)
01749     for(yi=0; yi<img->dimy; yi++) for(xi=0; xi<img->dimx; xi++)
01750       *fptr++=img->m[zi][yi][xi][frame_index];
01751 
01752   /* Write subheader and data, and set the rest of subheader contents */
01753   fptr=fdata; ret=-1;
01754   for(zi=0; zi<img->dimz; zi++, fptr+=img->dimx*img->dimy) {
01755     /* Create new matrix id (i.e. matnum) */
01756     matrixId=mat_numcod(frame_to_write, img->planeNumber[zi], 1, 0, 0);
01757     if(img->type==IMG_TYPE_RAW) {
01758       scan_header.frame_start_time=
01759         (int)temp_roundf(1000.*img->start[frame_index]);
01760       scan_header.frame_duration=
01761         (int)temp_roundf(1000.*(img->end[frame_index]-img->start[frame_index]));
01762       scan_header.prompts=temp_roundf(img->prompts[frame_index]);
01763       scan_header.delayed=temp_roundf(img->randoms[frame_index]);
01764       ret=ecat63WriteScanMatrix(fp, matrixId, &scan_header, fptr);
01765     } else {
01766       image_header.frame_start_time=
01767         (int)temp_roundf(1000.*img->start[frame_index]);
01768       image_header.frame_duration=
01769         (int)temp_roundf(1000.*(img->end[frame_index]-img->start[frame_index]));
01770       if(img->decayCorrected)
01771         image_header.decay_corr_fctr=img->decayCorrFactor[frame_index];
01772       ret=ecat63WriteImageMatrix(fp, matrixId, &image_header, fptr);
01773     }
01774     if(ret!=0) break;
01775   } /* next plane*/
01776   free(fdata); fclose(fp);
01777   if(ret) return STATUS_DISKFULL;
01778 
01779   return STATUS_OK;
01780 }
01781 /*****************************************************************************/
01782 
01783 /*****************************************************************************/
01790 void imgSetEcat63SHeader(IMG *img, void *h) {
01791   ECAT63_imageheader *image_header;
01792   ECAT63_scanheader *scan_header;
01793 
01794   if(img->type==IMG_TYPE_RAW) {
01795     scan_header=(ECAT63_scanheader*)h;
01796     scan_header->data_type=VAX_I2;
01797     scan_header->dimension_1=img->dimx;
01798     scan_header->dimension_2=img->dimy;
01799     scan_header->frame_duration_sec=1;
01800     scan_header->scale_factor=1.0;
01801     scan_header->frame_start_time=0;
01802     scan_header->frame_duration=1000;
01803     scan_header->loss_correction_fctr=1.0;
01804   } else {
01805     image_header=(ECAT63_imageheader*)h;
01806     image_header->data_type=VAX_I2;
01807     image_header->num_dimensions=2;
01808     image_header->dimension_1=img->dimx;
01809     image_header->dimension_2=img->dimy;
01810     image_header->recon_scale=img->zoom;
01811     image_header->quant_scale=1.0;
01812     image_header->slice_width=img->sizez/10.;
01813     image_header->pixel_size=img->sizex/10.;
01814     image_header->frame_start_time=0;
01815     image_header->frame_duration=1000;
01816     image_header->plane_eff_corr_fctr=1.0;
01817     image_header->decay_corr_fctr=1.0;
01818     image_header->loss_corr_fctr=1.0;
01819     image_header->ecat_calibration_fctr=1.0;
01820     image_header->well_counter_cal_fctr=1.0;
01821     image_header->quant_units=imgUnitToEcat6(img);
01822   }
01823 }
01824 /*****************************************************************************/
01825 
01826 /*****************************************************************************/