img_e7.c
Go to the documentation of this file.
00001 /******************************************************************************
00002 
00003   Copyright (c) 2007,2010 Turku PET Centre
00004 
00005   Library:     img_e7.c
00006   Description: ECAT 7 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-02-27 Vesa Oikonen
00027     Functions moved from imgfile.c.
00028     Bug correction in imgWrite2DEcat7():
00029     prompts were set after the header was written.
00030   2007-03-25 VO
00031     Added functions imgReadEcat7Header(), imgEcat7Supported(),
00032     imgReadEcat7Frame(), imgReadEcat7FirstFrame(), imgGetEcat7Fileformat(),
00033     imgWriteEcat7Frame(), and imgSetEcat7SHeader().
00034   2007-04-03 VO
00035     Added support for ECAT 7 polar maps.
00036   2007-17-07 Harri Merisaari
00037     Modified timezone correction in for ANSI compatibility and for 
00038     other timezones in imgGetEcat7MHeader
00039   2007-09-07 VO
00040     The rest of strcpy's replaced by strncpy's in imgGetEcat7MHeader(),
00041     most importantly study_description, where overflow caused accidental 
00042     change of IMG type. 
00043   2007-09-10 VO
00044     Return value of localtime() is checked.
00045   2007-09-12 VO
00046     Bug correction in imgReadEcat7Frame(): did not set 2D image plane numbers.
00047   2010-02-12 VO
00048     Timezone correction added also into imgSetEcat7MHeader().
00049 
00050 
00051 ******************************************************************************/
00052 #include <stdio.h>
00053 #include <stdlib.h>
00054 #include <unistd.h>
00055 #include <math.h>
00056 #include <string.h>
00057 #include <time.h>
00058 /*****************************************************************************/
00059 #include "petc99.h"
00060 #include "swap.h"
00061 #include "halflife.h"
00062 /*****************************************************************************/
00063 #include "include/img.h"
00064 #include "include/ecat7.h"
00065 #include "include/imgmax.h"
00066 #include "include/imgdecay.h"
00067 #include "include/imgfile.h"
00068 /*****************************************************************************/
00069 
00070 /*****************************************************************************/
00083 int imgReadEcat7(const char *fname, IMG *img) {
00084   FILE *fp;
00085   int ret, m, i, fi, pi, xi, yi, frame, plane, prev_frame, prev_plane;
00086   int dimx, dimy, dimz, planeNr, frameNr, blkNr=0, pxlNr;
00087   ECAT7_mainheader main_header;
00088   ECAT7_imageheader image_header;
00089   ECAT7_2Dscanheader scan2d_header;
00090   ECAT7_scanheader scan_header;
00091   ECAT7_polmapheader polmap_header;
00092   ECAT7_MATRIXLIST mlist;
00093   ECAT7_Matval matval;
00094   float *fdata=NULL, *fptr;
00095 
00096 
00097   if(IMG_TEST) printf("imgReadEcat7(%s, *img)\n", fname);
00098   /* Check the arguments */
00099   if(fname==NULL) return(1);
00100   if(img==NULL || img->status!=IMG_STATUS_INITIALIZED) {
00101     imgSetStatus(img, STATUS_FAULT); return(2);}
00102 
00103   /* Open file for read */
00104   if((fp=fopen(fname, "rb"))==NULL) {imgSetStatus(img, STATUS_NOFILE); return(3);}
00105   
00106   /* Read main header */
00107   ret=ecat7ReadMainheader(fp, &main_header);
00108   if(ret) {fclose(fp); imgSetStatus(img, STATUS_UNKNOWNFORMAT); return(4);}
00109   /* Check for magic number */
00110   if(strncmp(main_header.magic_number, ECAT7V_MAGICNR, 7)!=0) {
00111     fclose(fp); imgSetStatus(img, STATUS_UNKNOWNFORMAT); return(4);
00112   }
00113   
00114   /* Check if file_type is supported */
00115   if(imgEcat7Supported(&main_header)==0) {
00116     fclose(fp); imgSetStatus(img, STATUS_UNSUPPORTED); return(5);
00117   }
00118 
00119   /* Read matrix list */
00120   ecat7InitMatlist(&mlist);
00121   ret=ecat7ReadMatlist(fp, &mlist);
00122   if(ret || mlist.matrixNr<1 || ecat7CheckMatlist(&mlist)) {
00123     fclose(fp); imgSetStatus(img, STATUS_INVALIDMATLIST); return(6);}
00124   ecat7SortMatlistByPlane(&mlist);
00125   if(IMG_TEST>2) ecat7PrintMatlist(&mlist);
00126 
00127   /* Calculate the number of planes, frames and gates */
00128   /* Check that all planes have equal nr of frames (gates) */
00129   /* and check that frames (gates) are consequentally numbered */
00130   prev_plane=plane=-1; prev_frame=frame=-1; frameNr=planeNr=0; ret=0;
00131   for(m=0; m<mlist.matrixNr; m++) {
00132     ecat7_id_to_val(mlist.matdir[m].id, &matval); plane=matval.plane;
00133     if(main_header.num_frames>=main_header.num_gates) frame=matval.frame;
00134     else frame=matval.gate;
00135     if(plane!=prev_plane) {frameNr=1; planeNr++;}
00136     else {frameNr++; if(prev_frame>0 && frame!=prev_frame+1) {ret=1; break;}}
00137     prev_plane=plane; prev_frame=frame;
00138     /* Calculate and check the size of one data matrix */
00139     if(m==0) {
00140       blkNr=mlist.matdir[m].endblk-mlist.matdir[m].strtblk;
00141     } else if(blkNr!=mlist.matdir[m].endblk-mlist.matdir[m].strtblk) {
00142       ret=2; break;
00143     }
00144   } /* next matrix */
00145   if(IMG_TEST>2) printf("frameNr=%d planeNr=%d\n", frameNr, planeNr);
00146   if(ret==1 || (frameNr*planeNr != mlist.matrixNr)) {
00147     fclose(fp); imgSetStatus(img, STATUS_MISSINGMATRIX); ecat7EmptyMatlist(&mlist); return(7);}
00148   if(ret==2) {
00149     fclose(fp); imgSetStatus(img, STATUS_VARMATSIZE); ecat7EmptyMatlist(&mlist); return(8);}
00150 
00151   /* Read the first subheader to get planeNr from volumes and to get x&y dim */
00152   m=0; dimz=1; imgSetStatus(img, STATUS_NOSUBHEADER);
00153   switch(main_header.file_type) {
00154     case ECAT7_IMAGE8:
00155     case ECAT7_IMAGE16:
00156     case ECAT7_VOLUME8:
00157     case ECAT7_VOLUME16:
00158       img->type=IMG_TYPE_IMAGE;
00159       ret=ecat7ReadImageheader(fp, mlist.matdir[m].strtblk, &image_header);
00160       dimx=image_header.x_dimension; dimy=image_header.y_dimension;
00161       if(image_header.num_dimensions>2 && image_header.z_dimension>1)
00162         planeNr=dimz=image_header.z_dimension;
00163       break;
00164     case ECAT7_2DSCAN:
00165       img->type=IMG_TYPE_RAW;
00166       ret=ecat7Read2DScanheader(fp, mlist.matdir[m].strtblk, &scan2d_header);
00167       dimx=scan2d_header.num_r_elements; dimy=scan2d_header.num_angles;
00168       if(scan2d_header.num_dimensions>2 && scan2d_header.num_z_elements>1)
00169         planeNr=dimz=scan2d_header.num_z_elements;
00170       break;
00171     case ECAT7_3DSCAN:
00172     case ECAT7_3DSCAN8:
00173     case ECAT7_3DSCANFIT:
00174       img->type=IMG_TYPE_RAW;
00175       ret=ecat7ReadScanheader(fp, mlist.matdir[m].strtblk, &scan_header);
00176       dimx=scan_header.num_r_elements; dimy=scan_header.num_angles;
00177       for(i=dimz=0; i<64; i++) dimz+=scan_header.num_z_elements[i];
00178       planeNr=dimz;
00179       /*if(scan_header.axial_compression!=0) {img->statmsg=imgmsg[STATUS_UNSUPPORTEDAXIALCOMP]; ret=-1;}*/
00180       break;
00181     case ECAT7_POLARMAP:
00182       img->type=IMG_TYPE_POLARMAP;
00183       ret=ecat7ReadPolmapheader(fp, mlist.matdir[m].strtblk, &polmap_header);
00184       planeNr=dimz=dimy=1; dimx=0;
00185       for(i=0; i<polmap_header.num_rings; i++)
00186         dimx+=polmap_header.sectors_per_ring[i];
00187       break;
00188     default: dimx=dimy=dimz=planeNr=0; ret=-1;
00189   }
00190   pxlNr=dimx*dimy;
00191   if(ret || pxlNr<1 || planeNr<1) {
00192     fclose(fp);  ecat7EmptyMatlist(&mlist); return(9);}
00193   imgSetStatus(img, STATUS_OK);
00194 
00195   /* Allocate memory for IMG data */
00196   ret=imgAllocate(img, planeNr, dimy, dimx, frameNr);
00197   if(ret) {
00198     fclose(fp); imgSetStatus(img, STATUS_NOMEMORY);
00199     ecat7EmptyMatlist(&mlist); return(11);
00200   }
00201   /* Copy information from mainheader */
00202   imgGetEcat7MHeader(img, &main_header);
00203   /* Set fileFormat */
00204   switch(main_header.file_type) {
00205     case ECAT7_IMAGE8:
00206     case ECAT7_IMAGE16:
00207       img->_fileFormat=IMG_E7_2D; break;
00208     case ECAT7_VOLUME8:
00209     case ECAT7_VOLUME16:
00210       img->_fileFormat=IMG_E7; break;
00211     case ECAT7_2DSCAN:
00212       img->_fileFormat=IMG_E7_2D; break;
00213     case ECAT7_3DSCAN:
00214     case ECAT7_3DSCAN8:
00215     case ECAT7_3DSCANFIT:
00216       img->_fileFormat=IMG_E7; break;
00217     case ECAT7_POLARMAP:
00218       img->_fileFormat=IMG_POLARMAP; break;
00219     default:
00220       img->_fileFormat=IMG_UNKNOWN; break;
00221   }
00222 
00223   if(dimz>1) {
00224     /* Read ECAT volume matrices */
00225     fi=0;
00226     for(m=0; m<mlist.matrixNr; m++) {
00227       /* Get matrix values */
00228       ecat7_id_to_val(mlist.matdir[m].id, &matval);
00229       /* Read subheader and data */
00230       if(img->type==IMG_TYPE_IMAGE)
00231         ret=ecat7ReadImageMatrix(fp, mlist.matdir[m].strtblk,
00232               mlist.matdir[m].endblk, &image_header, &fdata);
00233       else
00234         ret=ecat7ReadScanMatrix(fp, mlist.matdir[m].strtblk,
00235               mlist.matdir[m].endblk, &scan_header, &fdata);
00236       if(ret || fdata==NULL) {
00237         if(IMG_TEST) printf("ecat7ReadXMatrix()=%d\n%s\n", ret, ecat7errmsg);
00238         fclose(fp); imgSetStatus(img, STATUS_NOMATRIX); ecat7EmptyMatlist(&mlist); return(13);}
00239       /* Copy subheader information */
00240       if(img->type==IMG_TYPE_POLARMAP) {
00241         img->_dataType=polmap_header.data_type;
00242         img->start[fi]=polmap_header.frame_start_time/1000.;
00243         img->end[fi]=img->start[fi]+polmap_header.frame_duration/1000.;
00244         img->mid[fi]=0.5*(img->start[fi]+img->end[fi]);
00245         img->sizex=0.001*polmap_header.pixel_size;
00246       } else if(img->type==IMG_TYPE_IMAGE) {
00247         img->_dataType=image_header.data_type;
00248         img->start[fi]=image_header.frame_start_time/1000.;
00249         img->end[fi]=img->start[fi]+image_header.frame_duration/1000.;
00250         img->mid[fi]=0.5*(img->start[fi]+img->end[fi]);
00251         if(image_header.decay_corr_fctr>1.0) {
00252           img->decayCorrFactor[fi]=image_header.decay_corr_fctr;
00253           img->decayCorrected=1;
00254         }
00255         img->zoom=image_header.recon_zoom;
00256         img->sizex=10.*image_header.x_pixel_size;
00257         img->sizey=10.*image_header.y_pixel_size;
00258         img->sizez=10.*image_header.z_pixel_size;
00259         img->resolutionx=10.*image_header.x_resolution;
00260         img->resolutiony=10.*image_header.y_resolution;
00261         img->resolutionz=10.*image_header.z_resolution;
00262       } else {
00263         img->_dataType=scan_header.data_type;
00264         img->start[fi]=scan_header.frame_start_time/1000.;
00265         img->end[fi]=img->start[fi]+scan2d_header.frame_duration/1000.;
00266         img->mid[fi]=0.5*(img->start[fi]+img->end[fi]);
00267         if(scan_header.x_resolution>0.0)
00268           img->sampleDistance=10.0*scan_header.x_resolution;
00269         else
00270           img->sampleDistance=10.0*main_header.bin_size;
00271         /* also, correct for dead-time */
00272         if(scan_header.deadtime_correction_factor>0.0)
00273           for(i=0, fptr=fdata; i<dimz*pxlNr; i++, fptr++)
00274             *fptr*=scan_header.deadtime_correction_factor;
00275         img->prompts[fi]=(float)scan_header.prompts;
00276         img->randoms[fi]=scan_header.delayed;
00277       }
00278       /* Copy matrix data through volume planes */
00279       for(pi=0; pi<dimz; pi++) {
00280         for(yi=0, fptr=fdata+pi*pxlNr; yi<dimy; yi++) for(xi=0; xi<dimx; xi++)
00281           img->m[pi][yi][xi][fi]= *fptr++;
00282       }
00283       free(fdata); fi++;
00284     } /* next matrix */
00285     /* Set plane numbers */
00286     for(pi=0; pi<dimz; pi++) img->planeNumber[pi]=pi+1;
00287   } else {
00288     /* Read separate matrices */
00289     prev_plane=plane=-1; prev_frame=frame=-1; pi=fi=-1;
00290     for(m=0; m<mlist.matrixNr; m++) {
00291       ecat7_id_to_val(mlist.matdir[m].id, &matval); plane=matval.plane;
00292       if(main_header.num_frames>=main_header.num_gates) frame=matval.frame;
00293       else frame=matval.gate;
00294       if(plane!=prev_plane) {fi=0; pi++;} else fi++;
00295       /* Read subheader and data */
00296       if(img->type==IMG_TYPE_POLARMAP)
00297         ret=ecat7ReadPolarmapMatrix(fp, mlist.matdir[m].strtblk,
00298               mlist.matdir[m].endblk, &polmap_header, &fdata);
00299       else if(img->type==IMG_TYPE_IMAGE)
00300         ret=ecat7ReadImageMatrix(fp, mlist.matdir[m].strtblk,
00301               mlist.matdir[m].endblk, &image_header, &fdata);
00302       else
00303         ret=ecat7Read2DScanMatrix(fp, mlist.matdir[m].strtblk,
00304               mlist.matdir[m].endblk, &scan2d_header, &fdata);
00305       if(ret || fdata==NULL) {
00306         fclose(fp); imgSetStatus(img, STATUS_NOMATRIX); ecat7EmptyMatlist(&mlist); return(13);}
00307       /* Copy subheader information */
00308       if(fi==0) img->planeNumber[pi]=plane;
00309       if(img->type==IMG_TYPE_POLARMAP) {
00310         img->_dataType=polmap_header.data_type;
00311         img->start[fi]=polmap_header.frame_start_time/1000.;
00312         img->end[fi]=img->start[fi]+polmap_header.frame_duration/1000.;
00313         img->mid[fi]=0.5*(img->start[fi]+img->end[fi]);
00314         img->sizex=0.001*polmap_header.pixel_size;
00315       } else if(img->type==IMG_TYPE_IMAGE) {
00316         img->_dataType=image_header.data_type;
00317         img->start[fi]=image_header.frame_start_time/1000.;
00318         img->end[fi]=img->start[fi]+image_header.frame_duration/1000.;
00319         img->mid[fi]=0.5*(img->start[fi]+img->end[fi]);
00320         if(image_header.decay_corr_fctr>1.0) {
00321           img->decayCorrFactor[fi]=image_header.decay_corr_fctr;
00322           img->decayCorrected=1;
00323         }
00324         img->zoom=image_header.recon_zoom;
00325         img->sizex=10.*image_header.x_pixel_size;
00326         img->sizey=10.*image_header.y_pixel_size;
00327         img->sizez=10.*image_header.z_pixel_size;
00328         img->resolutionx=10.*image_header.x_resolution;
00329         img->resolutiony=10.*image_header.y_resolution;
00330         img->resolutionz=10.*image_header.z_resolution;
00331       } else {
00332         img->_dataType=scan2d_header.data_type;
00333         img->start[fi]=scan2d_header.frame_start_time/1000.;
00334         img->end[fi]=img->start[fi]+scan2d_header.frame_duration/1000.;
00335         img->mid[fi]=0.5*(img->start[fi]+img->end[fi]);
00336         if(scan_header.x_resolution>0.0)
00337           img->sampleDistance=10.0*scan_header.x_resolution;
00338         else
00339           img->sampleDistance=10.0*main_header.bin_size;
00340         /* also, correct for dead-time */
00341         if(scan2d_header.deadtime_correction_factor>0.0)
00342           for(i=0, fptr=fdata; i<pxlNr; i++, fptr++)
00343             *fptr*=scan2d_header.deadtime_correction_factor;
00344         img->prompts[fi]=(float)scan_header.prompts;
00345         img->randoms[fi]=scan_header.delayed;
00346       }
00347       /* Copy matrix data */
00348       for(yi=0, fptr=fdata; yi<dimy; yi++) for(xi=0; xi<dimx; xi++)
00349         img->m[pi][yi][xi][fi]= *fptr++;
00350       free(fdata);
00351       /* prepare for the next matrix */
00352       prev_plane=plane; prev_frame=frame;
00353     } /* next matrix */
00354   }
00355   fclose(fp); ecat7EmptyMatlist(&mlist);
00356 
00357   /* Calibrate */
00358   if(main_header.ecat_calibration_factor>0.0)
00359     for(pi=0; pi<img->dimz; pi++)
00360       for(yi=0; yi<img->dimy; yi++) for(xi=0; xi<img->dimx; xi++)
00361         for(fi=0; fi<img->dimt; fi++)
00362           img->m[pi][yi][xi][fi]*=main_header.ecat_calibration_factor;
00363 
00364   imgSetStatus(img, STATUS_OK);
00365   return(0);
00366 }
00367 /*****************************************************************************/
00368 
00369 /*****************************************************************************/
00380 int imgWriteEcat7(const char *fname, IMG *img) {
00381   ECAT7_mainheader main_header;
00382   ECAT7_imageheader image_header;
00383   ECAT7_scanheader scan_header;
00384   FILE *fp;
00385   int fi, pi, xi, yi, pxlNr, matrixId, ret;
00386   float *fdata, *fptr;
00387 
00388 
00389   if(IMG_TEST) printf("imgWriteEcat7(%s, *img)\n", fname);
00390   if(IMG_TEST>1 && ECAT7_TEST==0) ECAT7_TEST=1;
00391   /* Check the arguments */
00392   if(fname==NULL) return(1);
00393   if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) {
00394     imgSetStatus(img, STATUS_FAULT); return(2);}
00395   if(img->type!=IMG_TYPE_RAW && img->type!=IMG_TYPE_IMAGE) {
00396     imgSetStatus(img, STATUS_FAULT); return(2);}
00397 
00398   /* Initiate headers */
00399   memset(&main_header, 0, sizeof(ECAT7_mainheader));
00400   memset(&image_header,0, sizeof(ECAT7_imageheader));
00401   memset(&scan_header, 0, sizeof(ECAT7_scanheader));
00402 
00403   /* Set main header */
00404   imgSetEcat7MHeader(img, &main_header);
00405   main_header.bin_size=img->sampleDistance/10.0;
00406 
00407   /* Allocate memory for matrix float data */
00408   pxlNr=img->dimx*img->dimy*img->dimz;
00409   fdata=(float*)malloc(pxlNr*sizeof(float));
00410   if(fdata==NULL) {imgSetStatus(img, STATUS_NOMEMORY); return(3);}
00411 
00412   /* Open file, write main header and initiate matrix list */
00413   fp=ecat7Create(fname, &main_header);
00414   if(fp==NULL) {free(fdata); imgSetStatus(img, STATUS_NOWRITEPERM); return(6);}
00415 
00416   /* Set (most of) subheader contents */
00417   if(img->type==IMG_TYPE_RAW) {
00418     scan_header.x_resolution=img->sampleDistance/10.0;
00419     scan_header.num_dimensions=4;
00420     if(img->dimz==239) {
00421       scan_header.num_z_elements[0]=63;
00422       scan_header.num_z_elements[1]=106;
00423       scan_header.num_z_elements[2]=70;
00424     } else {
00425       scan_header.num_z_elements[0]=img->dimz;
00426     }
00427     scan_header.storage_order=1;
00428     scan_header.data_type=ECAT7_SUNI2;
00429     scan_header.num_r_elements=img->dimx;
00430     scan_header.num_angles=img->dimy;
00431   } else if(img->type==IMG_TYPE_IMAGE) {
00432     image_header.num_dimensions=3;
00433     image_header.z_dimension=img->dimz;
00434     image_header.data_type=ECAT7_SUNI2;
00435     image_header.x_dimension=img->dimx;
00436     image_header.y_dimension=img->dimy;
00437     image_header.recon_zoom=img->zoom;
00438     image_header.x_pixel_size=0.1*img->sizex;
00439     image_header.y_pixel_size=0.1*img->sizey;
00440     image_header.z_pixel_size=0.1*img->sizez;
00441     image_header.x_resolution=0.1*img->resolutionx;
00442     image_header.y_resolution=0.1*img->resolutiony;
00443     image_header.z_resolution=0.1*img->resolutionz;
00444   }
00445 
00446   /* Write each matrix */
00447   for(fi=0; fi<img->dimt; fi++) {
00448 
00449     /* Create new matrix id (i.e. matnum) */
00450     matrixId=ecat7_val_to_id(fi+1, 1, 1, 0, 0);
00451 
00452     /* Copy matrix pixel values to fdata */
00453     fptr=fdata;
00454     for(pi=0; pi<img->dimz; pi++)
00455       for(yi=0; yi<img->dimy; yi++) for(xi=0; xi<img->dimx; xi++)
00456         *fptr++=img->m[pi][yi][xi][fi];
00457 
00458     /* Write subheader and data */
00459     fptr=fdata;
00460     if(img->type==IMG_TYPE_RAW) {
00461       scan_header.frame_start_time=(int)temp_roundf(1000.*img->start[fi]);
00462       scan_header.frame_duration=(int)temp_roundf(1000.*(img->end[fi]-img->start[fi]));
00463       scan_header.prompts=temp_roundf(img->prompts[fi]);
00464       scan_header.delayed=temp_roundf(img->randoms[fi]);
00465       /*ecat7PrintScanheader(&scan_header, stdout);*/
00466       ret=ecat7WriteScanMatrix(fp, matrixId, &scan_header, fptr);
00467     } else if(img->type==IMG_TYPE_IMAGE) {
00468       image_header.frame_start_time=(int)temp_roundf(1000.*img->start[fi]);
00469       image_header.frame_duration=(int)temp_roundf(1000.*(img->end[fi]-img->start[fi]));
00470       image_header.decay_corr_fctr=img->decayCorrFactor[fi];
00471       /*ecat7PrintImageheader(&image_header, stdout);*/
00472       ret=ecat7WriteImageMatrix(fp, matrixId, &image_header, fptr);
00473     } else {free(fdata); fclose(fp); imgSetStatus(img, STATUS_UNSUPPORTED); return(8);}
00474     if(ret) {
00475       if(IMG_TEST) {printf("matrixId=%d ret=%d\n", matrixId, ret);}
00476       free(fdata); fclose(fp); imgSetStatus(img, STATUS_DISKFULL); return(7);
00477     }
00478 
00479   } /* next matrix */
00480   free(fdata); fclose(fp);
00481 
00482   imgSetStatus(img, STATUS_OK);
00483   return(0);
00484 }
00485 /*****************************************************************************/
00486 
00487 /*****************************************************************************/
00498 int imgWrite2DEcat7(const char *fname, IMG *img) {
00499   ECAT7_mainheader main_header;
00500   ECAT7_imageheader image_header;
00501   ECAT7_2Dscanheader scan2d_header;
00502   FILE *fp;
00503   int fi, pi, xi, yi, pxlNr, matrixId, ret;
00504   float *fdata, *fptr;
00505 
00506 
00507   if(IMG_TEST) printf("imgWrite2DEcat7(%s, *img)\n", fname);
00508   if(IMG_TEST>1 && ECAT7_TEST==0) ECAT7_TEST=1;
00509   /* Check the arguments */
00510   if(fname==NULL) {return(1);}
00511   if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) {
00512     imgSetStatus(img, STATUS_FAULT); return(2);}
00513 
00514   /* Initiate headers */
00515   memset(&main_header, 0, sizeof(ECAT7_mainheader));
00516   memset(&image_header,0, sizeof(ECAT7_imageheader));
00517   memset(&scan2d_header, 0, sizeof(ECAT7_2Dscanheader));
00518 
00519   /* Set main header */
00520   imgSetEcat7MHeader(img, &main_header);
00521   main_header.bin_size=img->sampleDistance/10.0;
00522   if(img->type==IMG_TYPE_RAW) main_header.file_type=ECAT7_2DSCAN;
00523   else main_header.file_type=ECAT7_IMAGE16;
00524   main_header.num_planes=img->dimz;
00525 
00526   /* Allocate memory for matrix float data */
00527   pxlNr=img->dimx*img->dimy;
00528   fdata=(float*)malloc(pxlNr*sizeof(float));
00529   if(fdata==NULL) {imgSetStatus(img, STATUS_NOMEMORY); return(3);}
00530 
00531   /* Open file, write main header and initiate matrix list */
00532   fp=ecat7Create(fname, &main_header);
00533   if(fp==NULL) {free(fdata); imgSetStatus(img, STATUS_NOWRITEPERM); return(6);}
00534 
00535   /* Set (most of) subheader contents */
00536   if(img->type==IMG_TYPE_RAW) {
00537     scan2d_header.num_dimensions=2;
00538     scan2d_header.num_z_elements=1;
00539     scan2d_header.data_type=ECAT7_SUNI2;
00540     scan2d_header.num_r_elements=img->dimx;
00541     scan2d_header.num_angles=img->dimy;
00542   } else if(img->type==IMG_TYPE_IMAGE) {
00543     image_header.num_dimensions=2;
00544     image_header.z_dimension=1;
00545     image_header.data_type=ECAT7_SUNI2;
00546     image_header.x_dimension=img->dimx;
00547     image_header.y_dimension=img->dimy;
00548     image_header.recon_zoom=img->zoom;
00549     image_header.x_pixel_size=0.1*img->sizex;
00550     image_header.y_pixel_size=0.1*img->sizey;
00551     image_header.z_pixel_size=0.1*img->sizez;
00552     image_header.x_resolution=0.1*img->resolutionx;
00553     image_header.y_resolution=0.1*img->resolutiony;
00554     image_header.z_resolution=0.1*img->resolutionz;
00555   }
00556 
00557   /* Write each matrix */
00558   for(fi=0; fi<img->dimt; fi++) for(pi=0; pi<img->dimz; pi++) {
00559 
00560     /* Create new matrix id (i.e. matnum) */
00561     matrixId=ecat7_val_to_id(fi+1, img->planeNumber[pi], 1, 0, 0);
00562 
00563     /* Copy matrix pixel values to fdata */
00564     fptr=fdata;
00565     for(yi=0; yi<img->dimy; yi++) for(xi=0; xi<img->dimx; xi++)
00566       *fptr++=img->m[pi][yi][xi][fi];
00567 
00568     /* Write subheader and data */
00569     fptr=fdata;
00570     if(img->type==IMG_TYPE_RAW) {
00571       scan2d_header.frame_start_time=(int)temp_roundf(1000.*img->start[fi]);
00572       scan2d_header.frame_duration=(int)temp_roundf(1000.*(img->end[fi]-img->start[fi]));
00573       scan2d_header.prompts=temp_roundf(img->prompts[fi]);
00574       scan2d_header.delayed=temp_roundf(img->randoms[fi]);
00575       ret=ecat7Write2DScanMatrix(fp, matrixId, &scan2d_header, fptr);
00576     } else if(img->type==IMG_TYPE_IMAGE) {
00577       image_header.frame_start_time=(int)temp_roundf(1000.*img->start[fi]);
00578       image_header.frame_duration=(int)temp_roundf(1000.*(img->end[fi]-img->start[fi]));
00579       image_header.decay_corr_fctr=img->decayCorrFactor[fi];
00580       ret=ecat7WriteImageMatrix(fp, matrixId, &image_header, fptr);
00581     } else {free(fdata); fclose(fp); imgSetStatus(img, STATUS_UNSUPPORTED); return(8);}
00582     if(ret) {
00583       if(IMG_TEST) {printf("matrixId=%d ret=%d\n", matrixId, ret);}
00584       free(fdata); fclose(fp); imgSetStatus(img, STATUS_DISKFULL); return(7);
00585     }
00586 
00587   } /* next matrix */
00588   free(fdata); fclose(fp);
00589 
00590   imgSetStatus(img, STATUS_OK);
00591   return(0);
00592 }
00593 /*****************************************************************************/
00594 
00595 /*****************************************************************************/
00606 int imgWritePolarmap(const char *fname, IMG *img) {
00607   ECAT7_mainheader main_header;
00608   ECAT7_polmapheader polmap_header;
00609   FILE *fp;
00610   int fi, pi, xi, yi, pxlNr, matrixId, ret;
00611   float *fdata, *fptr;
00612 
00613 
00614   if(IMG_TEST) printf("imgWritePolarmap(%s, *img)\n", fname);
00615   if(IMG_TEST>1 && ECAT7_TEST==0) ECAT7_TEST=1;
00616   /* Check the arguments */
00617   if(fname==NULL) return(1);
00618   if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) {
00619     imgSetStatus(img, STATUS_FAULT); return(2);}
00620   if(img->type!=IMG_TYPE_POLARMAP) {
00621     imgSetStatus(img, STATUS_FAULT); return(2);}
00622 
00623   /* Initiate headers */
00624   memset(&main_header, 0, sizeof(ECAT7_mainheader));
00625   memset(&polmap_header,0, sizeof(ECAT7_polmapheader));
00626 
00627   /* Set main header */
00628   imgSetEcat7MHeader(img, &main_header);
00629   main_header.bin_size=img->sampleDistance/10.0;
00630 
00631   /* Allocate memory for matrix float data */
00632   pxlNr=img->dimx*img->dimy*img->dimz;
00633   fdata=(float*)malloc(pxlNr*sizeof(float));
00634   if(fdata==NULL) {imgSetStatus(img, STATUS_NOMEMORY); return(3);}
00635 
00636   /* Open file, write main header and initiate matrix list */
00637   fp=ecat7Create(fname, &main_header);
00638   if(fp==NULL) {free(fdata); imgSetStatus(img, STATUS_NOWRITEPERM); return(6);}
00639 
00640   /* Set (most of) subheader contents */
00641   imgSetEcat7SHeader(img, &polmap_header);
00642 
00643   /* Write each matrix */
00644   for(fi=0; fi<img->dimt; fi++) {
00645 
00646     /* Create new matrix id (i.e. matnum) */
00647     matrixId=ecat7_val_to_id(fi+1, 1, 1, 0, 0);
00648 
00649     /* Copy matrix pixel values to fdata */
00650     fptr=fdata;
00651     for(pi=0; pi<img->dimz; pi++)
00652       for(yi=0; yi<img->dimy; yi++) for(xi=0; xi<img->dimx; xi++)
00653         *fptr++=img->m[pi][yi][xi][fi];
00654 
00655     /* Write subheader and data */
00656     fptr=fdata;
00657     polmap_header.frame_start_time=(int)temp_roundf(1000.*img->start[fi]);
00658     polmap_header.frame_duration=(int)temp_roundf(1000.*(img->end[fi]-img->start[fi]));
00659     /*ecat7PrintImageheader(&image_header, stdout);*/
00660     ret=ecat7WritePolarmapMatrix(fp, matrixId, &polmap_header, fptr);
00661     if(ret) {
00662       if(IMG_TEST) {printf("matrixId=%d ret=%d\n", matrixId, ret);}
00663       free(fdata); fclose(fp); imgSetStatus(img, STATUS_DISKFULL); return(7);
00664     }
00665 
00666   } /* next matrix */
00667   free(fdata); fclose(fp);
00668 
00669   imgSetStatus(img, STATUS_OK);
00670   return(0);
00671 }
00672 /*****************************************************************************/
00673 
00674 /*****************************************************************************/
00681 void imgGetEcat7MHeader(IMG *img, ECAT7_mainheader *h) {
00682 #ifndef __STRICT_ANSI__
00683   struct tm *scanStartTime;
00684   time_t lt;
00685 #endif
00686   
00687   img->scanner=h->system_type;
00688   imgUnitFromEcat7(img, h);
00689   strncpy(img->radiopharmaceutical, h->radiopharmaceutical, 32);
00690   img->isotopeHalflife=h->isotope_halflife;
00691   img->scanStart=h->scan_start_time;
00692 #ifndef __STRICT_ANSI__
00693   /* refresh timezone information (might be out of date) */
00694   tzset();
00695   /* correct scan start time in hours */
00696   lt = (time_t)img->scanStart;
00697   scanStartTime = localtime(&lt);
00698   if(scanStartTime!=NULL) {
00699     scanStartTime->tm_hour -= (timezone/3600);
00700     img->scanStart = mktime(scanStartTime);
00701   }
00702 #endif
00703   img->axialFOV=10.0*h->distance_scanned;
00704   img->transaxialFOV=10.0*h->transaxial_fov;
00705   strncpy(img->studyNr, h->study_type, MAX_STUDYNR_LEN);
00706   strncpy(img->patientName, h->patient_name, 32);
00707   strncpy(img->patientID, h->patient_id, 16);
00708   img->sizez=10.0*h->plane_separation;
00709   switch(h->file_type) {
00710     case ECAT7_IMAGE8:
00711     case ECAT7_IMAGE16:
00712     case ECAT7_VOLUME8:
00713     case ECAT7_VOLUME16: img->type=IMG_TYPE_IMAGE; break;
00714     case ECAT7_POLARMAP: img->type=IMG_TYPE_POLARMAP; break;
00715     default: img->type=IMG_TYPE_RAW;
00716   }
00717   img->orientation=h->patient_orientation;
00718   strncpy(img->studyDescription, h->study_description, 32);
00719   strncpy(img->userProcessCode, h->user_process_code, 10);
00720   img->userProcessCode[10]=(char)0;
00721   /* If valid study number is found in user_process_code, then take it */  
00722   if(!img->studyNr[0] && studynr_validity_check(img->userProcessCode))
00723     strcpy(img->studyNr, img->userProcessCode);
00724 }
00725 /*****************************************************************************/
00726 
00727 /*****************************************************************************/
00734 void imgSetEcat7MHeader(IMG *img, ECAT7_mainheader *h) {
00735 #ifndef __STRICT_ANSI__
00736   struct tm *scanStartTime;
00737   time_t lt;
00738 #endif
00739 
00740   h->sw_version=72;
00741   if(img->type==IMG_TYPE_POLARMAP) {
00742     strcpy(h->magic_number, ECAT7V_MAGICNR);
00743     h->file_type=ECAT7_POLARMAP;
00744   } else if(img->type==IMG_TYPE_RAW) {
00745     strcpy(h->magic_number, ECAT7S_MAGICNR);
00746     if(img->_fileFormat==IMG_E7_2D) h->file_type=ECAT7_2DSCAN;
00747     else h->file_type=ECAT7_3DSCAN;
00748   } else {
00749     strcpy(h->magic_number, ECAT7V_MAGICNR);
00750     if(img->_fileFormat==IMG_E7_2D) h->file_type=ECAT7_IMAGE16;
00751     else h->file_type=ECAT7_VOLUME16;
00752   }
00753   h->system_type=img->scanner;
00754   h->scan_start_time=img->scanStart;
00755 #ifndef __STRICT_ANSI__
00756   /* refresh timezone information (might be out of date) */
00757   tzset();
00758   /* correct scan start time in hours */
00759   lt = (time_t)h->scan_start_time;
00760   scanStartTime = localtime(&lt);
00761   if(scanStartTime!=NULL) {
00762     scanStartTime->tm_hour += (timezone/3600);
00763     h->scan_start_time = mktime(scanStartTime);
00764   }
00765 #endif
00766   h->isotope_halflife=img->isotopeHalflife;
00767   imgUnitToEcat7(img, h);
00768   h->ecat_calibration_factor=1.0;
00769   h->transaxial_fov=img->transaxialFOV/10.0;
00770   h->num_planes=img->dimz; /* h->num_planes=1; */
00771   h->num_frames=img->dimt;  
00772   h->num_gates=1;
00773   h->num_bed_pos=0;
00774   h->distance_scanned=img->axialFOV/10.0;
00775   h->plane_separation=img->sizez/10.0;
00776   strncpy(h->radiopharmaceutical, img->radiopharmaceutical, 32);
00777   strcpy(h->isotope_name, imgIsotope(img));
00778   strcpy(h->study_type, img->studyNr);
00779   strcpy(h->patient_name, img->patientName);
00780   strcpy(h->patient_id, img->patientID);
00781   h->patient_orientation=img->orientation;
00782   strcpy(h->study_description, img->studyDescription);
00783   strncpy(h->user_process_code, img->userProcessCode, 10);
00784 }
00785 /*****************************************************************************/
00786 
00787 /*****************************************************************************/
00794 int imgGetEcat7Fileformat(ECAT7_mainheader *h) {
00795   int fileFormat=IMG_UNKNOWN;
00796   switch(h->file_type) {
00797     case ECAT7_IMAGE8:
00798     case ECAT7_IMAGE16:
00799       fileFormat=IMG_E7_2D; break;
00800     case ECAT7_VOLUME8:
00801     case ECAT7_VOLUME16:
00802       fileFormat=IMG_E7; break;
00803     case ECAT7_2DSCAN:
00804       fileFormat=IMG_E7_2D; break;
00805     case ECAT7_3DSCAN:
00806     case ECAT7_3DSCAN8:
00807     case ECAT7_3DSCANFIT:
00808       fileFormat=IMG_E7; break;
00809     case ECAT7_POLARMAP:
00810       fileFormat=IMG_POLARMAP; break;
00811     default:
00812       fileFormat=IMG_UNKNOWN; break;
00813   }
00814   return fileFormat;
00815 }
00816 /*****************************************************************************/
00827 int imgReadEcat7Header(const char *fname, IMG *img) {
00828   FILE *fp;
00829   int ret, m, i;
00830   int planeNr, frameNr, blkNr=0;
00831   ECAT7_mainheader main_header;
00832   ECAT7_imageheader image_header;
00833   ECAT7_2Dscanheader scan2d_header;
00834   ECAT7_scanheader scan_header;
00835   ECAT7_polmapheader polmap_header;
00836   ECAT7_MATRIXLIST mlist;
00837 
00838 
00839   if(IMG_TEST) printf("\nimgReadEcat7Header(%s, *img)\n", fname);
00840   
00841   /* Check the arguments */
00842   if(img==NULL) return STATUS_FAULT;
00843   if(img->status!=IMG_STATUS_INITIALIZED) return STATUS_FAULT;
00844   imgSetStatus(img, STATUS_FAULT);
00845   if(fname==NULL) return STATUS_FAULT;
00846     
00847   /* Open the file */
00848   if((fp=fopen(fname, "rb")) == NULL) return STATUS_NOFILE;
00849 
00850   /* Read main header */
00851   ret=ecat7ReadMainheader(fp, &main_header);
00852   if(ret) {fclose(fp); return STATUS_NOMAINHEADER;}
00853   /* Check for magic number */
00854   if(strncmp(main_header.magic_number, ECAT7V_MAGICNR, 7)!=0) {
00855     fclose(fp); return STATUS_UNKNOWNFORMAT;}
00856   /* Check if file_type is supported */
00857   if(imgEcat7Supported(&main_header)==0) {fclose(fp); return STATUS_UNSUPPORTED;}
00858   /* Copy main header information into IMG; sets also img.type */
00859   imgGetEcat7MHeader(img, &main_header);
00860   if(IMG_TEST>7) printf("img.type := %d\n", img->type);
00861   img->_fileFormat=imgGetEcat7Fileformat(&main_header);
00862   if(IMG_TEST>7) printf("img._fileFormat := %d\n", img->_fileFormat);
00863   if(img->_fileFormat==IMG_UNKNOWN) {fclose(fp); return STATUS_UNSUPPORTED;}
00864 
00865   /* Read matrix list */
00866   ecat7InitMatlist(&mlist);
00867   ret=ecat7ReadMatlist(fp, &mlist);
00868   if(ret || mlist.matrixNr<1 || ecat7CheckMatlist(&mlist)) {
00869     fclose(fp); return STATUS_INVALIDMATLIST;}
00870   /* Make sure that plane and frame numbers are continuous */
00871   ecat7GatherMatlist(&mlist, 1, 1, 1, 1);
00872   /* Get plane and frame numbers and check that volume is full */
00873   ret=ecat7GetPlaneAndFrameNr(&mlist, &main_header, &planeNr, &frameNr);
00874   if(ret) {ecat7EmptyMatlist(&mlist); fclose(fp); return ret;}
00875   img->dimz=planeNr;
00876   img->dimt=frameNr;
00877   /* Get and check the size of data matrices */
00878   ret=ecat7GetMatrixBlockSize(&mlist, &blkNr);
00879   if(ret) {ecat7EmptyMatlist(&mlist); fclose(fp); return ret;}
00880 
00881   /* Read one subheader */
00882   if(IMG_TEST>5) printf("main_header.file_type := %d\n", main_header.file_type);
00883   m=0;
00884   switch(main_header.file_type) {
00885     case ECAT7_IMAGE8:
00886     case ECAT7_IMAGE16:
00887     case ECAT7_VOLUME8:
00888     case ECAT7_VOLUME16:
00889       ret=ecat7ReadImageheader(fp, mlist.matdir[m].strtblk, &image_header);
00890       break;
00891     case ECAT7_2DSCAN:
00892       ret=ecat7Read2DScanheader(fp, mlist.matdir[m].strtblk, &scan2d_header);
00893       break;
00894     case ECAT7_3DSCAN:
00895     case ECAT7_3DSCAN8:
00896     case ECAT7_3DSCANFIT:
00897       ret=ecat7ReadScanheader(fp, mlist.matdir[m].strtblk, &scan_header);
00898       break;
00899     case ECAT7_POLARMAP:
00900       ret=ecat7ReadPolmapheader(fp, mlist.matdir[m].strtblk, &polmap_header);
00901       break;
00902     default: ret=-1;
00903   }
00904   /* Free locally allocated memory and close the file */
00905   ecat7EmptyMatlist(&mlist); fclose(fp);
00906   /* Check whether subheader was read */
00907   if(ret) return STATUS_NOSUBHEADER;
00908 
00909   /* Get the following information in the subheader:
00910      dimensions x, y and z; datatype;
00911      image decay correction on/off, zoom, pixel size and resolution;
00912      sinogram sample distance.
00913    */
00914   switch(main_header.file_type) {
00915     case ECAT7_IMAGE8:
00916     case ECAT7_IMAGE16:
00917     case ECAT7_VOLUME8:
00918     case ECAT7_VOLUME16:
00919       img->dimx=image_header.x_dimension; img->dimy=image_header.y_dimension;
00920       if(image_header.num_dimensions>2 && image_header.z_dimension>1)
00921         /*planeNr=*/ img->dimz=image_header.z_dimension;
00922       img->_dataType=image_header.data_type;
00923       if(image_header.decay_corr_fctr>1.0) img->decayCorrected=1;
00924       img->zoom=image_header.recon_zoom;
00925       img->sizex=10.*image_header.x_pixel_size;
00926       img->sizey=10.*image_header.y_pixel_size;
00927       img->sizez=10.*image_header.z_pixel_size;
00928       img->resolutionx=10.*image_header.x_resolution;
00929       img->resolutiony=10.*image_header.y_resolution;
00930       img->resolutionz=10.*image_header.z_resolution;
00931       break;
00932     case ECAT7_2DSCAN:
00933       img->dimx=scan2d_header.num_r_elements; img->dimy=scan2d_header.num_angles;
00934       if(scan2d_header.num_dimensions>2 && scan2d_header.num_z_elements>1)
00935         planeNr=img->dimz=scan2d_header.num_z_elements;
00936       img->_dataType=scan2d_header.data_type;
00937       if(scan2d_header.x_resolution>0.0)
00938         img->sampleDistance=10.0*scan2d_header.x_resolution;
00939       else
00940         img->sampleDistance=10.0*main_header.bin_size;
00941       break;
00942     case ECAT7_3DSCAN:
00943     case ECAT7_3DSCAN8:
00944     case ECAT7_3DSCANFIT:
00945       img->dimx=scan_header.num_r_elements; img->dimy=scan_header.num_angles;
00946       for(i=img->dimz=0; i<64; i++) img->dimz+=scan_header.num_z_elements[i];
00947       /* planeNr=img->dimz; */
00948       img->_dataType=scan_header.data_type;
00949       if(scan_header.x_resolution>0.0)
00950         img->sampleDistance=10.0*scan_header.x_resolution;
00951       else
00952         img->sampleDistance=10.0*main_header.bin_size;
00953       break;
00954     case ECAT7_POLARMAP:
00955       img->dimy=img->dimz=1;
00956       img->polarmap_num_rings=polmap_header.num_rings;
00957       if(img->polarmap_num_rings>MAX_POLARMAP_NUM_RINGS)
00958         return STATUS_INVALIDPOLARMAP;
00959       for(i=0; i<img->polarmap_num_rings; i++) {
00960         img->polarmap_sectors_per_ring[i]=polmap_header.sectors_per_ring[i];
00961         img->polarmap_ring_position[i]=polmap_header.ring_position[i];
00962         img->polarmap_ring_angle[i]=polmap_header.ring_angle[i];
00963       }
00964       img->polarmap_start_angle=polmap_header.start_angle;
00965       for(i=0, img->dimx=0; i<img->polarmap_num_rings; i++)
00966         img->dimx+=img->polarmap_sectors_per_ring[i];
00967       /* pixel_size actually contains volume, in cubic cm */
00968       img->sizex=img->sizey=img->sizez=0.001*polmap_header.pixel_size;
00969       break;
00970   }
00971 
00972   imgSetStatus(img, STATUS_OK);
00973   return STATUS_OK;
00974 }
00975 /*****************************************************************************/
00976 
00977 /*****************************************************************************/
00984 int imgEcat7Supported(ECAT7_mainheader *h) {
00985   if(h==NULL) return(0);
00986   if(h->file_type==ECAT7_VOLUME8) return(1);
00987   if(h->file_type==ECAT7_VOLUME16) return(1);
00988   if(h->file_type==ECAT7_IMAGE8) return(1);
00989   if(h->file_type==ECAT7_IMAGE16) return(1);
00990   if(h->file_type==ECAT7_2DSCAN) return(1);
00991   if(h->file_type==ECAT7_3DSCAN) return(1);
00992   if(h->file_type==ECAT7_3DSCAN8) return(1);
00993   if(h->file_type==ECAT7_3DSCANFIT) return(1);
00994   if(h->file_type==ECAT7_POLARMAP) return(1);
00995   return(0);
00996 }
00997 /*****************************************************************************/
00998 
00999 /*****************************************************************************/
01008 int imgReadEcat7FirstFrame(const char *fname, IMG *img) {
01009   int ret=0;
01010 
01011   if(IMG_TEST) printf("\nimgReadEcat7FirstFrame(%s, *img)\n", fname);
01012   /* Check the input */
01013   if(img==NULL) return STATUS_FAULT;
01014   if(img->status!=IMG_STATUS_INITIALIZED) return STATUS_FAULT;
01015   imgSetStatus(img, STATUS_FAULT);
01016   if(fname==NULL) return STATUS_FAULT;
01017 
01018   /* Read header information from file */
01019   ret=imgReadEcat7Header(fname, img); if(ret) return(ret);
01020   if(IMG_TEST>3) imgInfo(img);
01021 
01022   /* Allocate memory for one frame */
01023   img->dimt=1;
01024   ret=imgAllocate(img, img->dimz, img->dimy, img->dimx, img->dimt);
01025   if(ret) return STATUS_NOMEMORY;
01026 
01027   /* Read the first frame */
01028   ret=imgReadEcat7Frame(fname, 1, img, 0); if(ret) return(ret); 
01029 
01030   /* All went well */
01031   imgSetStatus(img, STATUS_OK);
01032   return STATUS_OK;
01033 }
01034 /*****************************************************************************/
01035 
01036 /*****************************************************************************/
01051 int imgReadEcat7Frame(const char *fname, int frame_to_read, IMG *img, int frame_index) {
01052   FILE *fp;
01053   int ret, m, i, pi, xi, yi, frame, plane, seqplane, pxlNr;
01054   /* int blkNr=0; */
01055   ECAT7_mainheader main_header;
01056   ECAT7_imageheader image_header;
01057   ECAT7_2Dscanheader scan2d_header;
01058   ECAT7_scanheader scan_header;
01059   ECAT7_polmapheader polmap_header;
01060   ECAT7_MATRIXLIST mlist;
01061   ECAT7_Matval matval;
01062   float *fdata=NULL, *fptr;
01063 
01064 
01065   if(IMG_TEST) printf("\nimgReadEcat7Frame(%s, %d, *img, %d)\n",
01066     fname, frame_to_read, frame_index);
01067     
01068   /* Check the input */
01069   if(img==NULL) return STATUS_FAULT;
01070   if(img->status!=IMG_STATUS_OCCUPIED) return STATUS_FAULT;
01071   if(fname==NULL) return STATUS_FAULT;
01072   if(frame_index<0 || frame_index>img->dimt-1) return STATUS_FAULT;
01073   if(frame_to_read<1) return STATUS_FAULT;
01074 
01075   /* Open file for read */
01076   if((fp=fopen(fname, "rb")) == NULL) return STATUS_NOFILE;
01077 
01078   /* Read main header */
01079   ret=ecat7ReadMainheader(fp, &main_header);
01080   if(ret) {fclose(fp); return STATUS_NOMAINHEADER;}
01081   /* Check for magic number */
01082   if(strncmp(main_header.magic_number, ECAT7V_MAGICNR, 7)!=0) {
01083     fclose(fp); return STATUS_UNKNOWNFORMAT;}
01084   /* Check if file_type is supported */
01085   if(imgEcat7Supported(&main_header)==0) {
01086     fclose(fp); return STATUS_UNSUPPORTED;}
01087   /* Read matrix list and nr and sort it */
01088   ecat7InitMatlist(&mlist);
01089   ret=ecat7ReadMatlist(fp, &mlist);
01090   if(ret) {fclose(fp); return STATUS_NOMATLIST;}
01091   if(mlist.matrixNr<=0 || ecat7CheckMatlist(&mlist)) {
01092     fclose(fp); ecat7EmptyMatlist(&mlist); return STATUS_INVALIDMATLIST;}
01093   /* Make sure that plane and frame numbers are continuous */
01094   ecat7GatherMatlist(&mlist, 1, 1, 1, 1);
01095   ecat7SortMatlistByFrame(&mlist); /* printf("matlist sorted\n"); */
01096   /* Calculate and check the size of one data matrix */
01097   /*ret=ecat7GetMatrixBlockSize(&mlist, &blkNr); if(ret) {fclose(fp); return ret;}*/
01098 
01099   /* Read all matrices that belong to the required frame */
01100   /*blkNr=-1;*/
01101   ret=0; seqplane=-1; pxlNr=img->dimx*img->dimy;
01102   for(m=0; m<mlist.matrixNr; m++) {
01103     /* get frame and plane */
01104     ecat7_id_to_val(mlist.matdir[m].id, &matval); plane=matval.plane;
01105     if(main_header.num_frames>=main_header.num_gates) frame=matval.frame;
01106     else frame=matval.gate; /* printf("frame=%d plane=%d\n", frame, plane); */
01107     if(frame!=frame_to_read) continue; /* do not process other frames */
01108     /*if(img->dimz>1) seqplane++; else seqplane=plane-1; */
01109     if(img->_fileFormat==IMG_E7_2D) seqplane=plane-1; else seqplane++;
01110     
01111     /* Read subheader and data */
01112     if(IMG_TEST>4) printf("reading matrix %d,%d\n", frame, plane);
01113     if(img->type==IMG_TYPE_IMAGE) { /* 2D or 3D image */
01114       ret=ecat7ReadImageMatrix(fp, mlist.matdir[m].strtblk,
01115               mlist.matdir[m].endblk, &image_header, &fdata);
01116     } else if(img->type==IMG_TYPE_POLARMAP) {  /* polarmap */
01117       ret=ecat7ReadPolarmapMatrix(fp, mlist.matdir[m].strtblk,
01118               mlist.matdir[m].endblk, &polmap_header, &fdata);
01119     } else if(img->dimz>1) { /* 3D sinogram */
01120       ret=ecat7ReadScanMatrix(fp, mlist.matdir[m].strtblk,
01121             mlist.matdir[m].endblk, &scan_header, &fdata);
01122     } else { /* 2D sinogram */
01123       ret=ecat7Read2DScanMatrix(fp, mlist.matdir[m].strtblk,
01124             mlist.matdir[m].endblk, &scan2d_header, &fdata);
01125     }
01126     if(ret || fdata==NULL) {
01127       fclose(fp); ecat7EmptyMatlist(&mlist); return STATUS_NOMATRIX;}
01128 
01129     /* Copy information concerning this frame and make correction to data */
01130     if(img->type==IMG_TYPE_IMAGE) {
01131       img->start[frame_index]=image_header.frame_start_time/1000.;
01132       img->end[frame_index]=img->start[frame_index]+image_header.frame_duration/1000.;
01133       img->mid[frame_index]=0.5*(img->start[frame_index]+img->end[frame_index]);
01134       if(image_header.decay_corr_fctr>1.0) {
01135         img->decayCorrFactor[frame_index]=image_header.decay_corr_fctr;}
01136     } else if(img->type==IMG_TYPE_POLARMAP) {  /* polarmap */
01137       img->start[frame_index]=polmap_header.frame_start_time/1000.;
01138       img->end[frame_index]=img->start[frame_index]+polmap_header.frame_duration/1000.;
01139       img->mid[frame_index]=0.5*(img->start[frame_index]+img->end[frame_index]);
01140     } else if(img->dimz>1) {
01141       img->start[frame_index]=scan_header.frame_start_time/1000.;
01142       img->end[frame_index]=img->start[frame_index]+scan2d_header.frame_duration/1000.;
01143       img->mid[frame_index]=0.5*(img->start[frame_index]+img->end[frame_index]);
01144       /* also, correct for dead-time */
01145       if(scan_header.deadtime_correction_factor>0.0)
01146         for(i=0, fptr=fdata; i<img->dimz*pxlNr; i++, fptr++)
01147           *fptr*=scan_header.deadtime_correction_factor;
01148       img->prompts[frame_index]=(float)scan_header.prompts;
01149       img->randoms[frame_index]=scan_header.delayed;
01150     } else {
01151       img->start[frame_index]=scan2d_header.frame_start_time/1000.;
01152       img->end[frame_index]=img->start[frame_index]+scan2d_header.frame_duration/1000.;
01153       img->mid[frame_index]=0.5*(img->start[frame_index]+img->end[frame_index]);
01154       /* also, correct for dead-time */
01155       if(scan2d_header.deadtime_correction_factor>0.0)
01156         for(i=0, fptr=fdata; i<pxlNr; i++, fptr++)
01157           *fptr*=scan2d_header.deadtime_correction_factor;
01158       img->prompts[frame_index]=(float)scan2d_header.prompts;
01159       img->randoms[frame_index]=scan2d_header.delayed;
01160     }
01161     /* Copy matrix data through volume planes */
01162     if(img->_fileFormat!=IMG_E7_2D) {
01163     /* if(img->dimz>1) { */
01164       for(pi=0; pi<img->dimz; pi++) {
01165         if(IMG_TEST>5) printf("  putting data into m[%d][][][%d]\n", pi, frame_index);
01166         for(yi=0, fptr=fdata+pi*pxlNr; yi<img->dimy; yi++)
01167           for(xi=0; xi<img->dimx; xi++)
01168             img->m[pi][yi][xi][frame_index]= *fptr++;
01169       }
01170     } else {
01171         if(IMG_TEST>5) printf("  putting data into m[%d][][][%d]\n", seqplane, frame_index);
01172         for(yi=0, fptr=fdata; yi<img->dimy; yi++)
01173           for(xi=0; xi<img->dimx; xi++)
01174             img->m[seqplane][yi][xi][frame_index]= *fptr++;
01175         img->planeNumber[seqplane]=plane;
01176     }
01177     free(fdata);
01178     /* Calibrate */
01179     if(main_header.ecat_calibration_factor>0.0)
01180       for(pi=0; pi<img->dimz; pi++)
01181         for(yi=0; yi<img->dimy; yi++) for(xi=0; xi<img->dimx; xi++)
01182           img->m[pi][yi][xi][frame_index]*=main_header.ecat_calibration_factor;
01183   } /* next matrix */
01184   if(IMG_TEST>3) printf("end of matrices.\n");
01185   ecat7EmptyMatlist(&mlist);
01186   fclose(fp);
01187 
01188   /* seqplane is <0 only if this frame did not exist at all; return -1 */
01189   if(IMG_TEST>4) printf("last_seqplane := %d.\n", seqplane);
01190   if(seqplane<0) return STATUS_NOMATRIX;
01191 
01192   /* check that correct number of planes was read */
01193   if(seqplane>0 && (seqplane+1 != img->dimz)) return STATUS_MISSINGMATRIX;
01194 
01195   /* All went well */
01196   imgSetStatus(img, STATUS_OK);
01197   return STATUS_OK;
01198 }
01199 /*****************************************************************************/
01200 
01201 /*****************************************************************************/
01222 int imgWriteEcat7Frame(const char *fname, int frame_to_write, IMG *img, int frame_index) {
01223   IMG test_img;
01224   int ret=0, pxlNr, zi, xi, yi, matrixId;
01225   ECAT7_mainheader main_header;
01226   ECAT7_imageheader image_header;
01227   ECAT7_scanheader scan_header;
01228   ECAT7_2Dscanheader scan2d_header;
01229   ECAT7_polmapheader polmap_header;
01230   void *sub_header=NULL;
01231   FILE *fp;
01232   float *fdata=NULL, *fptr;
01233 
01234 
01235   if(IMG_TEST) printf("\nimgWriteEcat7Frame(%s, %d, *img, %d)\n",
01236     fname, frame_to_write, frame_index);
01237     
01238   /*
01239    *  Check the input 
01240    */
01241   if(fname==NULL) return STATUS_FAULT;
01242   if(img==NULL) return STATUS_FAULT;
01243   if(img->status!=IMG_STATUS_OCCUPIED) return STATUS_FAULT;
01244   if(frame_to_write<0) return STATUS_FAULT;
01245   if(frame_index<0 || frame_index>=img->dimt) return STATUS_FAULT;
01246   if(img->_fileFormat!=IMG_E7 &&
01247      img->_fileFormat!=IMG_POLARMAP &&
01248      img->_fileFormat!=IMG_E7_2D) return STATUS_FAULT;
01249 
01250   /* Initiate headers */
01251   memset(&main_header, 0, sizeof(ECAT7_mainheader));
01252   memset(&image_header,0, sizeof(ECAT7_imageheader));
01253   memset(&scan_header, 0, sizeof(ECAT7_scanheader));
01254   memset(&scan2d_header, 0, sizeof(ECAT7_2Dscanheader));
01255   memset(&polmap_header,0, sizeof(ECAT7_polmapheader));
01256   imgInit(&test_img);
01257 
01258 
01259   /*
01260    *  If file does not exist, then create it with new mainheader,
01261    *  and if it does exist, then read and check header information.
01262    *  Create or edit mainheader to contain correct frame nr.   
01263    *  In any case, leave file pointer open for write.   
01264    */
01265   if(access(fname, 0) == -1) { /* file does not exist */
01266 
01267     /* Set main header */
01268     imgSetEcat7MHeader(img, &main_header);
01269     main_header.bin_size=img->sampleDistance/10.0;
01270     if(frame_to_write==0) frame_to_write=1;
01271     main_header.num_frames=frame_to_write;
01272 
01273     /* Open file, write main header and initiate matrix list */
01274     fp=ecat7Create(fname, &main_header); if(fp==NULL) return STATUS_NOWRITEPERM;
01275 
01276   } else { /* file exists */
01277   
01278     /* Read header information for checking */
01279     ret=imgReadEcat7Header(fname, &test_img); if(ret!=0) return ret;
01280     /* Check that file format is the same */
01281     if(img->_fileFormat!=test_img._fileFormat || img->type!=test_img.type)
01282       return STATUS_WRONGFILETYPE;
01283     /* Check that matrix sizes are the same */
01284     if(img->dimz!=test_img.dimz || img->dimx!=test_img.dimx ||
01285        img->dimy!=test_img.dimy)
01286       return STATUS_VARMATSIZE;
01287     imgEmpty(&test_img);
01288 
01289     /* Open the file for read and write */
01290     if((fp=fopen(fname, "r+b"))==NULL) return STATUS_NOWRITEPERM;
01291 
01292     /* Read the mainheader, set new frame number, and write it back */
01293     if((ret=ecat7ReadMainheader(fp, &main_header))!=0) return STATUS_NOMAINHEADER;
01294     if(frame_to_write==0) frame_to_write=main_header.num_frames+1;
01295     if(main_header.num_frames<frame_to_write)
01296       main_header.num_frames=frame_to_write;
01297     if((ret=ecat7WriteMainheader(fp, &main_header))!=0) return STATUS_NOWRITEPERM;
01298     
01299   }
01300   if(IMG_TEST>2) {
01301     printf("frame_to_write := %d\n", frame_to_write);
01302   }
01303 
01304   /* Allocate memory for matrix float data */
01305   pxlNr=img->dimx*img->dimy*img->dimz; /* for 2D too! */
01306   fdata=(float*)malloc(pxlNr*sizeof(float));
01307   if(fdata==NULL) {fclose(fp); return STATUS_NOMEMORY;}
01308   
01309   /* Set matrix subheader */
01310   if(img->_fileFormat==IMG_E7) {
01311     if(img->type==IMG_TYPE_RAW) sub_header=(void*)&scan_header;
01312     else sub_header=&image_header;
01313   } else if(img->_fileFormat==IMG_E7_2D) {
01314     if(img->type==IMG_TYPE_RAW) sub_header=(void*)&scan_header;
01315     else sub_header=(void*)&image_header;
01316   } else if(img->_fileFormat==IMG_POLARMAP) {
01317     sub_header=(void*)&polmap_header;
01318   } else {fclose(fp); return STATUS_FAULT;}
01319   imgSetEcat7SHeader(img, sub_header);
01320 
01321   /* Copy matrix pixel values to fdata */
01322   fptr=fdata;
01323   for(zi=0; zi<img->dimz; zi++)
01324     for(yi=0; yi<img->dimy; yi++) for(xi=0; xi<img->dimx; xi++)
01325       *fptr++=img->m[zi][yi][xi][frame_index];
01326 
01327   /* Write subheader and data, and set the rest of subheader contents */
01328   fptr=fdata; ret=-1;
01329   if(img->_fileFormat==IMG_E7) {
01330     /* Create new matrix id (i.e. matnum) */
01331     matrixId=ecat7_val_to_id(frame_to_write, 1, 1, 0, 0);
01332     if(img->type==IMG_TYPE_RAW) {
01333       scan_header.frame_start_time=(int)temp_roundf(1000.*img->start[frame_index]);
01334       scan_header.frame_duration=
01335         (int)temp_roundf(1000.*(img->end[frame_index]-img->start[frame_index]));
01336       scan_header.prompts=temp_roundf(img->prompts[frame_index]);
01337       scan_header.delayed=temp_roundf(img->randoms[frame_index]);
01338       ret=ecat7WriteScanMatrix(fp, matrixId, &scan_header, fptr);
01339     } else {
01340       image_header.frame_start_time=(int)temp_roundf(1000.*img->start[frame_index]);
01341       image_header.frame_duration=
01342         (int)temp_roundf(1000.*(img->end[frame_index]-img->start[frame_index]));
01343       image_header.decay_corr_fctr=img->decayCorrFactor[frame_index];
01344       /*ecat7PrintImageheader(&image_header, stdout);*/
01345       ret=ecat7WriteImageMatrix(fp, matrixId, &image_header, fptr);
01346     }
01347   } else if(img->_fileFormat==IMG_E7_2D) {
01348     for(zi=0; zi<img->dimz; zi++, fptr+=img->dimx*img->dimy) {
01349       /* Create new matrix id (i.e. matnum) */
01350       matrixId=ecat7_val_to_id(frame_to_write, img->planeNumber[zi], 1, 0, 0);
01351       if(img->type==IMG_TYPE_RAW) {
01352         scan2d_header.frame_start_time=
01353     (int)temp_roundf(1000.*img->start[frame_index]);
01354         scan2d_header.frame_duration=
01355           (int)temp_roundf(1000.*(img->end[frame_index]-img->start[frame_index]));
01356         scan2d_header.prompts=temp_roundf(img->prompts[frame_index]);
01357         scan2d_header.delayed=temp_roundf(img->randoms[frame_index]);
01358         ret=ecat7Write2DScanMatrix(fp, matrixId, &scan2d_header, fptr);
01359       } else {
01360         image_header.frame_start_time=
01361     (int)temp_roundf(1000.*img->start[frame_index]);
01362         image_header.frame_duration=
01363     (int)temp_roundf(1000.*(img->end[frame_index]-img->start[frame_index]));
01364         image_header.decay_corr_fctr=img->decayCorrFactor[frame_index];
01365         ret=ecat7WriteImageMatrix(fp, matrixId, &image_header, fptr);
01366       }
01367       if(ret!=0) break;
01368     } /* next plane */
01369   } else if(img->_fileFormat==IMG_POLARMAP) {
01370     /* Create new matrix id (i.e. matnum) */
01371     matrixId=ecat7_val_to_id(frame_to_write, 1, 1, 0, 0);
01372     polmap_header.frame_start_time=(int)temp_roundf(1000.*img->start[frame_index]);
01373     polmap_header.frame_duration=
01374         (int)temp_roundf(1000.*(img->end[frame_index]-img->start[frame_index]));
01375     ret=ecat7WritePolarmapMatrix(fp, matrixId, &polmap_header, fptr);
01376   }
01377   free(fdata); fclose(fp);
01378   if(ret) return STATUS_DISKFULL;
01379 
01380   return STATUS_OK;
01381 }
01382 /*****************************************************************************/
01383 
01384 /*****************************************************************************/
01391 void imgSetEcat7SHeader(IMG *img, void *h) {
01392   ECAT7_imageheader *image_header;
01393   ECAT7_scanheader *scan_header;
01394   ECAT7_2Dscanheader *scan2d_header;
01395   ECAT7_polmapheader *polmap_header;
01396   int i;
01397 
01398   if(img->type==IMG_TYPE_POLARMAP) {
01399     polmap_header=(ECAT7_polmapheader*)h;
01400     polmap_header->data_type=ECAT7_SUNI2;
01401     polmap_header->num_rings=img->polarmap_num_rings;
01402     if(polmap_header->num_rings>32) polmap_header->num_rings=32;
01403     for(i=0; i<polmap_header->num_rings; i++) {
01404       polmap_header->sectors_per_ring[i]=img->polarmap_sectors_per_ring[i];
01405       polmap_header->ring_position[i]=img->polarmap_ring_position[i];
01406       polmap_header->ring_angle[i]=img->polarmap_ring_angle[i];
01407     }
01408     polmap_header->start_angle=258;
01409     polmap_header->pixel_size=1000.0*img->sizex;  
01410     polmap_header->quant_units=0; /* default, see main header */
01411   } else if(img->type==IMG_TYPE_RAW) {
01412     if(img->_fileFormat==IMG_E7_2D) {
01413       scan2d_header=(ECAT7_2Dscanheader*)h;
01414       scan2d_header->num_dimensions=2;
01415       scan2d_header->num_z_elements=1;
01416       scan2d_header->data_type=ECAT7_SUNI2;
01417       scan2d_header->num_r_elements=img->dimx;
01418       scan2d_header->num_angles=img->dimy;
01419     } else {
01420       scan_header=(ECAT7_scanheader*)h;
01421       scan_header->x_resolution=img->sampleDistance/10.0;
01422       scan_header->num_dimensions=4;
01423       if(img->dimz==239) {
01424         scan_header->num_z_elements[0]=63;
01425         scan_header->num_z_elements[1]=106;
01426         scan_header->num_z_elements[2]=70;
01427       } else {
01428         scan_header->num_z_elements[0]=img->dimz;
01429       }
01430       scan_header->storage_order=1;
01431       scan_header->data_type=ECAT7_SUNI2;
01432       scan_header->num_r_elements=img->dimx;
01433       scan_header->num_angles=img->dimy;
01434     }
01435   } else {
01436     image_header=(ECAT7_imageheader*)h;
01437     image_header->data_type=ECAT7_SUNI2;
01438     image_header->x_dimension=img->dimx;
01439     image_header->y_dimension=img->dimy;
01440     image_header->recon_zoom=img->zoom;
01441     image_header->x_pixel_size=0.1*img->sizex;
01442     image_header->y_pixel_size=0.1*img->sizey;
01443     image_header->z_pixel_size=0.1*img->sizez;
01444     image_header->x_resolution=0.1*img->resolutionx;
01445     image_header->y_resolution=0.1*img->resolutiony;
01446     image_header->z_resolution=0.1*img->resolutionz;
01447     if(img->_fileFormat==IMG_E7_2D) {
01448       image_header->num_dimensions=2;
01449       image_header->z_dimension=1;
01450     } else {
01451       image_header->num_dimensions=3;
01452       image_header->z_dimension=img->dimz;
01453     }
01454     for(i=0; i<49; i++) image_header->fill_user[i]=0;
01455   }
01456 }
01457 /*****************************************************************************/
01458 
01459 /*****************************************************************************/