My Project
img_e63.c
Go to the documentation of this file.
1 /******************************************************************************
2 
3  Copyright (c) 2007,2008 Turku PET Centre
4 
5  Library: img_e63.c
6  Description: ECAT 6.3 I/O routines for IMG data.
7 
8  This library is free software; you can redistribute it and/or
9  modify it under the terms of the GNU Lesser General Public
10  License as published by the Free Software Foundation; either
11  version 2.1 of the License, or (at your option) any later version.
12 
13  This library is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
16  See the GNU Lesser General Public License for more details:
17  http://www.gnu.org/copyleft/lesser.html
18 
19  You should have received a copy of the GNU Lesser General Public License
20  along with this library/program; if not, write to the Free Software
21  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 
23  Turku PET Centre, Turku, Finland, http://www.turkupetcentre.fi/
24 
25  Modification history:
26  2007-01-31 Vesa Oikonen
27  Functions moved from imgfile.c.
28  If valid study number is found in user_process_code, then take it from there.
29  Patient_id and study_description are read and written.
30  Prompts and randoms (delayed) are read and written.
31  2007-03-25 VO
32  Added functions imgGetEcat63MHeader(), imgSetEcat63MHeader(),
33  imgEcat63Supported(), imgGetEcat63Fileformat(), imgReadEcat63Header(),
34  imgReadEcat63FirstFrame(), imgReadEcat63Frame(), imgWriteEcat63Frame(),
35  and imgSetEcat63SHeader().
36  2007-09-10 VO
37  Return value of localtime() is checked.
38  2008-11-06 VO
39  study_description is copied with strncpy(), not strcpy().
40 
41 
42 
43 ******************************************************************************/
44 #include <stdio.h>
45 #include <stdlib.h>
46 #include <unistd.h>
47 #include <math.h>
48 #include <string.h>
49 #include <time.h>
50 /*****************************************************************************/
51 #include "petc99.h"
52 #include "swap.h"
53 #include "halflife.h"
54 /*****************************************************************************/
55 #include "include/img.h"
56 #include "include/ecat63.h"
57 #include "include/ecat7.h"
58 #include "include/imgmax.h"
59 #include "include/imgdecay.h"
60 #include "include/sif.h"
61 #include "include/imgfile.h"
62 /*****************************************************************************/
63 
64 /*****************************************************************************/
77 int ecat63ReadAllToImg(const char *fname, IMG *img) {
78  int i, j, m, ret, blkNr=0, dim_x=0, dim_y=0, pxlNr=0;
79  int frameNr, planeNr, del_nr=0;
80  int frame, plane, prev_frame, prev_plane, seqplane;
81  FILE *fp;
82  ECAT63_mainheader main_header;
83  MATRIXLIST mlist;
84  MatDir tmpMatdir;
85  Matval matval;
86  ECAT63_imageheader image_header;
87  ECAT63_scanheader scan_header;
88  ECAT63_attnheader attn_header;
89  ECAT63_normheader norm_header;
90  float scale=1.0;
91  short int *sptr;
92  char *mdata=NULL, *mptr;
93  int *iptr;
94  struct tm scanStart;
95 
96 
97  if(IMG_TEST) printf("ecat63ReadAllToImg(%s, *img)\n", fname);
98  /* Check the arguments */
99  if(fname==NULL || img==NULL || img->status!=IMG_STATUS_INITIALIZED) {
100  if(img) imgSetStatus(img, STATUS_FAULT);
101  return 1;
102  }
103 
104  /* Open input CTI file for read */
105  if((fp=fopen(fname, "rb")) == NULL) {
107  return 3;
108  }
109 
110  /* Read main header */
111  if((ret=ecat63ReadMainheader(fp, &main_header))) {
113  return 4;
114  }
115  if(IMG_TEST>5) ecat63PrintMainheader(&main_header, stdout);
116 
117  /* Read matrix list and nr */
118  ecat63InitMatlist(&mlist);
119  ret=ecat63ReadMatlist(fp, &mlist);
120  if(ret) {
122  return 5;
123  }
124  if(mlist.matrixNr<=0) {
125  strcpy(ecat63errmsg, "matrix list is empty");
126  return 5;
127  }
128  /* Sort matrix list */
129  for(i=0; i<mlist.matrixNr-1; i++) for(j=i+1; j<mlist.matrixNr; j++) {
130  if(mlist.matdir[i].matnum>mlist.matdir[j].matnum) {
131  tmpMatdir=mlist.matdir[i];
132  mlist.matdir[i]=mlist.matdir[j]; mlist.matdir[j]=tmpMatdir;
133  }
134  }
135  if(IMG_TEST>6) {
136  printf("matrix list after sorting:\n");
137  ecat63PrintMatlist(&mlist);
138  }
139 
140  /* Trim extra frames */
141  if(main_header.num_frames>0) {
142  del_nr=ecat63DeleteLateFrames(&mlist, main_header.num_frames);
143  if(IMG_TEST && del_nr>0)
144  printf(" %d entries in matrix list will not be used.\n", del_nr);
145  }
146  /* Calculate the number of planes, frames and gates */
147  /* Check that all planes have equal nr of frames (gates) */
148  /* and check that frames (gates) are consequentally numbered */
149  prev_plane=plane=-1; prev_frame=frame=-1; frameNr=planeNr=0; ret=0;
150  for(m=0; m<mlist.matrixNr; m++) if(mlist.matdir[m].matstat==1) {
151  /* get frame and plane */
152  mat_numdoc(mlist.matdir[m].matnum, &matval);
153  plane=matval.plane;
154  if(main_header.num_frames>=main_header.num_gates) frame=matval.frame;
155  else frame=matval.gate;
156  if(IMG_TEST>2) printf("frame=%d plane=%d\n", frame, plane);
157  /* did plane change? */
158  if(plane!=prev_plane) {
159  frameNr=1; planeNr++;
160  } else {
161  frameNr++;
162  /* In each plane, frame (gate) numbers must be continuous */
163  if(prev_frame>0 && frame!=prev_frame+1) {ret=1; break;}
164  }
165  prev_plane=plane; prev_frame=frame;
166  /* Calculate and check the size of one data matrix */
167  if(m==0) {
168  blkNr=mlist.matdir[m].endblk-mlist.matdir[m].strtblk;
169  } else if(blkNr!=mlist.matdir[m].endblk-mlist.matdir[m].strtblk) {
170  ret=2; break;
171  }
172  } /* next matrix */
173  if(IMG_TEST) printf("frameNr=%d planeNr=%d\n", frameNr, planeNr);
174  if(ret==1 || (frameNr*planeNr != mlist.matrixNr-del_nr)) {
175  strcpy(ecat63errmsg, "missing matrix");
176  ecat63EmptyMatlist(&mlist); fclose(fp);
177  return(6); /* this number is used in imgRead() */
178  }
179  if(ret==2) {
180  strcpy(ecat63errmsg, "matrix sizes are different");
181  ecat63EmptyMatlist(&mlist); fclose(fp); return(7);
182  }
183  /* Read x,y-dimensions from 1st matrix subheader */
184  m=0; if(main_header.file_type==IMAGE_DATA) {
185  ret=ecat63ReadImageheader(fp, mlist.matdir[m].strtblk, &image_header);
186  dim_x=image_header.dimension_1; dim_y=image_header.dimension_2;
187  } else if(main_header.file_type==RAW_DATA) {
188  ret=ecat63ReadScanheader(fp, mlist.matdir[m].strtblk, &scan_header);
189  dim_x=scan_header.dimension_1; dim_y=scan_header.dimension_2;
190  } else if(main_header.file_type==ATTN_DATA) {
191  ret=ecat63ReadAttnheader(fp, mlist.matdir[m].strtblk, &attn_header);
192  dim_x=attn_header.dimension_1; dim_y=attn_header.dimension_2;
193  } else if(main_header.file_type==NORM_DATA) {
194  ret=ecat63ReadNormheader(fp, mlist.matdir[m].strtblk, &norm_header);
195  dim_x=norm_header.dimension_1; dim_y=norm_header.dimension_2;
196  }
197  pxlNr=dim_x*dim_y;
198  if(ret) {
199  sprintf(ecat63errmsg, "cannot read matrix %u subheader", mlist.matdir[m].matnum);
200  ecat63EmptyMatlist(&mlist); fclose(fp); return(8);
201  }
202 
203  /* Allocate memory for ECAT data matrix */
204  if(IMG_TEST>1) printf("allocating memory for %d blocks\n", blkNr);
205  mdata=(char*)malloc(blkNr*MatBLKSIZE); if(mdata==NULL) {
206  strcpy(ecat63errmsg, "out of memory");
207  fclose(fp); ecat63EmptyMatlist(&mlist); return(8);
208  }
209  /* Allocate memory for whole img data */
210  ret=imgAllocate(img, planeNr, dim_y, dim_x, frameNr);
211  if(ret) {
212  sprintf(ecat63errmsg, "out of memory (%d)", ret);
213  fclose(fp); ecat63EmptyMatlist(&mlist); free(mdata); return(9);
214  }
215 
216  /* Fill img info with ECAT main and sub header information */
217  img->scanner=main_header.system_type;
218  img->unit=main_header.calibration_units; /* this continues below */
219  strncpy(img->radiopharmaceutical, main_header.radiopharmaceutical, 32);
220  img->isotopeHalflife=main_header.isotope_halflife;
221  memset(&scanStart, 0, sizeof(struct tm));
222  scanStart.tm_year=main_header.scan_start_year-1900;
223  scanStart.tm_mon=main_header.scan_start_month-1;
224  scanStart.tm_mday=main_header.scan_start_day; scanStart.tm_yday=0;
225  scanStart.tm_hour=main_header.scan_start_hour;
226  scanStart.tm_min=main_header.scan_start_minute;
227  scanStart.tm_sec=main_header.scan_start_second;
228  scanStart.tm_isdst=-1;
229  img->scanStart=mktime(&scanStart); if(img->scanStart==-1) img->scanStart=0;
230  img->axialFOV=10.*main_header.axial_fov;
231  img->transaxialFOV=10.*main_header.transaxial_fov;
232  strncpy(img->studyNr, main_header.study_name, MAX_STUDYNR_LEN);
233  strncpy(img->patientName, main_header.patient_name, 31);
234  strncpy(img->patientID, main_header.patient_id, 15);
235  img->sizez=10.*main_header.plane_separation;
236  if(main_header.file_type==IMAGE_DATA) {
237  img->type=IMG_TYPE_IMAGE;
238  if(img->unit<1) img->unit=image_header.quant_units;
239  img->zoom=image_header.recon_scale;
240  if(image_header.decay_corr_fctr>1.0) img->decayCorrected=1;
241  img->sizex=img->sizey=10.*image_header.pixel_size;
242  if(img->sizez<10.*image_header.slice_width)
243  img->sizez=10.*image_header.slice_width;
244  } else if(main_header.file_type==RAW_DATA) {
245  img->type=IMG_TYPE_RAW;
246  img->decayCorrected=0;
247  } else if(main_header.file_type==ATTN_DATA) {
248  img->type=IMG_TYPE_RAW;
249  img->decayCorrected=0;
250  } else if(main_header.file_type==NORM_DATA) {
251  img->type=IMG_TYPE_RAW;
252  img->decayCorrected=0;
253  }
254  strncpy(img->studyDescription, main_header.study_description, 32);
255  strncpy(img->userProcessCode, main_header.user_process_code, 10);
256  img->userProcessCode[10]=(char)0;
257  /* If valid study number is found in user_process_code, then take it */
258  if(!img->studyNr[0] && studynr_validity_check(img->userProcessCode))
259  strcpy(img->studyNr, img->userProcessCode);
260 
261  /* Set _fileFormat */
262  img->_fileFormat=IMG_E63;
263 
264  /* Read one ECAT matrix at a time and put them to img */
265  prev_plane=plane=-1; seqplane=-1;
266  for(m=0; m<mlist.matrixNr; m++) if(mlist.matdir[m].matstat==1) {
267  /* get plane */
268  mat_numdoc(mlist.matdir[m].matnum, &matval);
269  plane=matval.plane;
270  /* did plane change? */
271  if(plane!=prev_plane) {seqplane++; frame=1;} else frame++;
272  prev_plane=plane;
273  img->planeNumber[seqplane]=plane;
274  /* Read subheader */
275  if(main_header.file_type==IMAGE_DATA) {
276  ret=ecat63ReadImageheader(fp, mlist.matdir[m].strtblk, &image_header);
277  if(dim_x!=image_header.dimension_1 || dim_y!=image_header.dimension_2) ret=1;
278  scale=image_header.quant_scale;
279  if(image_header.ecat_calibration_fctr>0.0
280  && fabs(main_header.calibration_factor/image_header.ecat_calibration_fctr-1.0)>0.0001)
281  scale*=image_header.ecat_calibration_fctr;
282  if(img->unit==0 && image_header.quant_units>0) img->unit=image_header.quant_units;
283  img->_dataType=image_header.data_type;
284  img->start[frame-1]=image_header.frame_start_time/1000.;
285  img->end[frame-1]=img->start[frame-1]+image_header.frame_duration/1000.;
286  img->mid[frame-1]=0.5*(img->start[frame-1]+img->end[frame-1]);
287  if(image_header.decay_corr_fctr>1.0)
288  img->decayCorrFactor[frame-1]=image_header.decay_corr_fctr;
289  } else if(main_header.file_type==RAW_DATA) {
290  ret=ecat63ReadScanheader(fp, mlist.matdir[m].strtblk, &scan_header);
291  if(dim_x!=scan_header.dimension_1 || dim_y!=scan_header.dimension_2) ret=1;
292  scale=scan_header.scale_factor;
293  if(scan_header.loss_correction_fctr>0.0) scale*=scan_header.loss_correction_fctr;
294  img->_dataType=scan_header.data_type;
295  img->start[frame-1]=scan_header.frame_start_time/1000.;
296  img->end[frame-1]=img->start[frame-1]+scan_header.frame_duration/1000.;
297  img->mid[frame-1]=0.5*(img->start[frame-1]+img->end[frame-1]);
298  img->sampleDistance=10.0*scan_header.sample_distance;
299  img->prompts[frame-1]=(float)scan_header.prompts;
300  img->randoms[frame-1]=(float)scan_header.delayed;
301  } else if(main_header.file_type==ATTN_DATA) {
302  ret=ecat63ReadAttnheader(fp, mlist.matdir[m].strtblk, &attn_header);
303  if(dim_x!=attn_header.dimension_1 || dim_y!=attn_header.dimension_2) ret=1;
304  scale=attn_header.scale_factor;
305  img->_dataType=attn_header.data_type;
306  } else if(main_header.file_type==NORM_DATA) {
307  ret=ecat63ReadNormheader(fp, mlist.matdir[m].strtblk, &norm_header);
308  if(dim_x!=norm_header.dimension_1 || dim_y!=norm_header.dimension_2) ret=1;
309  scale=norm_header.scale_factor;
310  img->_dataType=norm_header.data_type;
311  } else img->_dataType=-1;
312  if(ret) {
313  sprintf(ecat63errmsg, "cannot read matrix %u subheader", mlist.matdir[m].matnum);
314  ecat63EmptyMatlist(&mlist); free(mdata); fclose(fp); return(10);
315  }
316  if(main_header.calibration_factor>0.0) scale*=main_header.calibration_factor;
317  if(IMG_TEST>2) printf("scale=%e datatype=%d\n", scale, img->_dataType);
318  /* Read the pixel data */
319  ret=ecat63ReadMatdata(fp, mlist.matdir[m].strtblk+1,
320  mlist.matdir[m].endblk-mlist.matdir[m].strtblk,
321  mdata, img->_dataType);
322  if(ret) {
323  strcpy(ecat63errmsg, "cannot read matrix data");
324  ecat63EmptyMatlist(&mlist); free(mdata); fclose(fp); return(11);
325  }
326  if(img->_dataType==BYTE_TYPE) {
327  for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr++) {
328  img->m[seqplane][i][j][frame-1]=scale*(float)(*mptr);
329  }
330  } else if(img->_dataType==VAX_I2 || img->_dataType==SUN_I2) {
331  for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr+=2) {
332  sptr=(short int*)mptr;
333  img->m[seqplane][i][j][frame-1]=scale*(float)(*sptr);
334  }
335  } else if(img->_dataType==VAX_I4 || img->_dataType==SUN_I4) {
336  for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr+=4) {
337  iptr=(int*)mptr;
338  img->m[seqplane][i][j][frame-1]=1.0; /*scale*(float)(*iptr);*/
339  }
340  } else if(img->_dataType==VAX_R4 || img->_dataType==IEEE_R4) {
341  for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr+=4) {
342  memcpy(&img->m[seqplane][i][j][frame-1], mptr, 4);
343  img->m[seqplane][i][j][frame-1]*=scale;
344  }
345  }
346  } /* next matrix */
347 
348  /* Set the unit */
349  imgUnitFromEcat(img, img->unit);
350 
351  /* Free memory and close file */
352  free(mdata);
353  ecat63EmptyMatlist(&mlist);
354  fclose(fp);
355 
356  /* For saving, only 2-byte VAX data types are allowed */
357  if(img->_dataType==VAX_I4 || img->_dataType==VAX_R4) img->_dataType=VAX_I2;
358 
359  return(0);
360 }
361 /*****************************************************************************/
362 
363 /*****************************************************************************/
374 int ecat63WriteAllImg(const char *fname, IMG *img) {
375  int frame, plane, n, i, j, m, ret=0;
376  float f, fmax, fmin, g, scale;
377  short int *sdata, *sptr, smin, smax;
378  FILE *fp;
379  ECAT63_mainheader main_header;
380  ECAT63_imageheader image_header;
381  ECAT63_scanheader scan_header;
382  struct tm *scanStart;
383 
384 
385  if(IMG_TEST) printf("ecat63WriteAllImg(%s, *img)\n", fname);
386  /* Check the arguments */
387  if(fname==NULL) {strcpy(ecat63errmsg, "invalid ECAT filename"); return(1);}
388  if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) {
389  strcpy(ecat63errmsg, "image data is empty"); return(2);}
390  if(img->_dataType<1) img->_dataType=VAX_I2;
391 
392  /* Initiate headers */
393  memset(&main_header, 0, sizeof(ECAT63_mainheader));
394  memset(&image_header,0, sizeof(ECAT63_imageheader));
395  memset(&scan_header, 0, sizeof(ECAT63_scanheader));
396 
397  /* Make a main header */
398  main_header.sw_version=2;
399  main_header.num_planes=img->dimz;
400  main_header.num_frames=img->dimt;
401  main_header.num_gates=1;
402  main_header.num_bed_pos=1;
403  if(img->type==IMG_TYPE_IMAGE) main_header.file_type=IMAGE_DATA;
404  else main_header.file_type=RAW_DATA;
405  main_header.data_type=img->_dataType;
406  if(img->scanner>0) main_header.system_type=img->scanner;
407  else main_header.system_type=ECAT63_SYSTEM_TYPE_DEFAULT;
408  main_header.calibration_factor=1.0;
409  main_header.axial_fov=img->axialFOV/10.0;
410  main_header.transaxial_fov=img->transaxialFOV/10.0;
411  main_header.plane_separation=img->sizez/10.0;
412  main_header.calibration_units=imgUnitToEcat6(img);
413  strncpy(main_header.radiopharmaceutical, img->radiopharmaceutical, 32);
414  scanStart=localtime(&img->scanStart);
415  if(scanStart!=NULL) {
416  main_header.scan_start_year=scanStart->tm_year+1900;
417  main_header.scan_start_month=scanStart->tm_mon+1;
418  main_header.scan_start_day=scanStart->tm_mday;
419  main_header.scan_start_hour=scanStart->tm_hour;
420  main_header.scan_start_minute=scanStart->tm_min;
421  main_header.scan_start_second=scanStart->tm_sec;
422  } else {
423  main_header.scan_start_year=1900;
424  main_header.scan_start_month=1;
425  main_header.scan_start_day=1;
426  main_header.scan_start_hour=0;
427  main_header.scan_start_minute=0;
428  main_header.scan_start_second=0;
429  }
430  main_header.isotope_halflife=img->isotopeHalflife;
431  strcpy(main_header.isotope_code, imgIsotope(img));
432  strcpy(main_header.study_name, img->studyNr);
433  strcpy(main_header.patient_name, img->patientName);
434  strcpy(main_header.patient_id, img->patientID);
435  strncpy(main_header.user_process_code, img->userProcessCode, 10);
436  strncpy(main_header.study_description, img->studyDescription, 32);
437  if(IMG_TEST) ecat63PrintMainheader(&main_header, stdout);
438 
439  /* Allocate memory for matrix data */
440  sdata=(short int*)malloc(img->dimx*img->dimy*sizeof(short int));
441  if(sdata==NULL) {strcpy(ecat63errmsg, "out of memory"); return(4);}
442 
443  /* Open output ECAT file */
444  fp=ecat63Create(fname, &main_header);
445  if(fp==NULL) {strcpy(ecat63errmsg, "cannot write ECAT file"); return(3);}
446 
447  /* Set the subheader contents, as far as possible */
448  switch(main_header.file_type) {
449  case RAW_DATA:
450  scan_header.data_type=main_header.data_type;
451  scan_header.dimension_1=img->dimx;
452  scan_header.dimension_2=img->dimy;
453  scan_header.frame_duration_sec=1;
454  scan_header.scale_factor=1.0;
455  scan_header.frame_start_time=0;
456  scan_header.frame_duration=1000;
457  scan_header.loss_correction_fctr=1.0;
458  /*if(IMG_TEST) ecat63PrintScanheader(&scan_header);*/
459  break;
460  case IMAGE_DATA:
461  image_header.data_type=main_header.data_type;
462  image_header.num_dimensions=2;
463  image_header.dimension_1=img->dimx;
464  image_header.dimension_2=img->dimy;
465  image_header.recon_scale=img->zoom;
466  image_header.quant_scale=1.0;
467  image_header.slice_width=img->sizez/10.;
468  image_header.pixel_size=img->sizex/10.;
469  image_header.frame_start_time=0;
470  image_header.frame_duration=1000;
471  image_header.plane_eff_corr_fctr=1.0;
472  image_header.decay_corr_fctr=1.0;
473  image_header.loss_corr_fctr=1.0;
474  image_header.ecat_calibration_fctr=1.0;
475  image_header.well_counter_cal_fctr=1.0;
476  image_header.quant_units=main_header.calibration_units;
477  /*if(IMG_TEST) ecat63PrintImageheader(&image_header);*/
478  break;
479  }
480 
481  /* Write one matrix at a time */
482  n=0;
483  for(plane=1; plane<=img->dimz; plane++) for(frame=1; frame<=img->dimt; frame++) {
484  /* Scale data to short ints */
485  /* Search min and max */
486  fmin=fmax=f=img->m[plane-1][0][0][frame-1];
487  for(i=0; i<img->dimy; i++) for(j=0; j<img->dimx; j++) {
488  f=img->m[plane-1][i][j][frame-1];
489  if(f>fmax) fmax=f; else if(f<fmin) fmin=f;
490  }
491  /* Calculate scaling factor */
492  if(fabs(fmin)>fabs(fmax)) g=fabs(fmin); else g=fabs(fmax);
493  if(g!=0) scale=32766./g; else scale=1.0;
494  /*printf("fmin=%e fmax=%e g=%e scale=%e\n", fmin, fmax, g, scale);*/
495  /* Scale matrix data to shorts */
496  sptr=sdata;
497  for(i=0; i<img->dimy; i++) for(j=0; j<img->dimx; j++) {
498  *sptr=(short int)temp_roundf(scale*img->m[plane-1][i][j][frame-1]);
499  sptr++;
500  }
501  /* Calculate and set subheader min&max and scale */
502  smin=(short int)temp_roundf(scale*fmin); smax=(short int)temp_roundf(scale*fmax);
503  if(main_header.file_type==RAW_DATA) {
504  scan_header.scan_min=smin; scan_header.scan_max=smax;
505  scan_header.scale_factor=1.0/scale;
506  scan_header.frame_start_time=(int)temp_roundf(1000.*img->start[frame-1]);
507  scan_header.frame_duration=
508  (int)temp_roundf(1000.*(img->end[frame-1]-img->start[frame-1]));
509  scan_header.sample_distance=(img->sampleDistance)/10.0;
510  scan_header.prompts=temp_roundf(img->prompts[frame-1]);
511  scan_header.delayed=temp_roundf(img->randoms[frame-1]);
512  } else if(main_header.file_type==IMAGE_DATA) {
513  image_header.image_min=smin; image_header.image_max=smax;
514  image_header.quant_scale=1.0/scale;
515  image_header.frame_start_time=(int)temp_roundf(1000.*img->start[frame-1]);
516  image_header.frame_duration=
517  (int)temp_roundf(1000.*(img->end[frame-1]-img->start[frame-1]));
518  /* Set decay correction factor */
519  if(img->decayCorrected)
520  image_header.decay_corr_fctr=img->decayCorrFactor[frame-1];
521  }
522  /* Write matrix data */
523  m=mat_numcod(frame, img->planeNumber[plane-1], 1, 0, 0);
524  /*m=mat_numcod(frame, plane, 1, 0, 0);*/
525  sptr=sdata;
526  if(IMG_TEST) printf(" writing matnum=%d\n", m);
527  if(main_header.file_type==RAW_DATA) {
528  if(IMG_TEST) ecat63PrintScanheader(&scan_header, stdout);
529  ret=ecat63WriteScan(fp, m, &scan_header, sptr);
530  } else if(main_header.file_type==IMAGE_DATA) {
531  if(IMG_TEST) ecat63PrintImageheader(&image_header, stdout);
532  ret=ecat63WriteImage(fp, m, &image_header, sptr);
533  }
534  if(ret) {
535  sprintf(ecat63errmsg, "cannot write data on pl%02d fr%02d (%d).",
536  plane, frame, ret);
537  fclose(fp); remove(fname); free(sdata);
538  return(9);
539  }
540  n++;
541  } /* next matrix */
542  if(IMG_TEST) printf(" %d matrices written in %s\n", n, fname);
543 
544  /* Close file and free memory */
545  fclose(fp); free(sdata);
546 
547  return(0);
548 }
549 /*****************************************************************************/
550 
551 /*****************************************************************************/
568 int ecat63ReadPlaneToImg(const char *fname, IMG *img) {
569  int i, j, m, ret, blkNr=0, dim_x=0, dim_y=0, pxlNr=0, del_nr=0;
570  int frameNr, planeNr, datatype=0, firstm=0, isFirst=0;
571  int frame, plane, prev_frame, prev_plane, prev_frameNr=0;
572  FILE *fp;
573  ECAT63_mainheader main_header;
574  MATRIXLIST mlist;
575  MatDir tmpMatdir;
576  Matval matval;
577  ECAT63_imageheader image_header;
578  ECAT63_scanheader scan_header;
579  ECAT63_attnheader attn_header;
580  ECAT63_normheader norm_header;
581  float scale=1.0;
582  short int *sptr;
583  char *mdata=NULL, *mptr;
584  int *iptr;
585  struct tm scanStart;
586 
587 
588  if(IMG_TEST) printf("ecat63ReadPlaneToImg(%s, *img)\n", fname);
589  /* Check the arguments */
590  if(fname==NULL) {strcpy(ecat63errmsg, "invalid ECAT filename"); return(2);}
591  if(img==NULL || img->status==IMG_STATUS_UNINITIALIZED) {
592  strcpy(ecat63errmsg, "image data not initialized"); return(2);}
593 
594  /* Open input CTI file for read */
595  if((fp=fopen(fname, "rb")) == NULL) {
596  strcpy(ecat63errmsg, "cannot open ECAT file"); return(3);}
597 
598  /* Read main header */
599  if((ret=ecat63ReadMainheader(fp, &main_header))) {
600  sprintf(ecat63errmsg, "cannot read main header (%d)", ret);
601  fclose(fp); return(4);
602  }
603  if(IMG_TEST) ecat63PrintMainheader(&main_header, stdout);
604 
605  /* Read matrix list and nr */
606  ecat63InitMatlist(&mlist);
607  ret=ecat63ReadMatlist(fp, &mlist);
608  if(ret) {
609  sprintf(ecat63errmsg, "cannot read matrix list (%d)", ret);
610  fclose(fp); return(5);
611  }
612  if(mlist.matrixNr<=0) {
613  strcpy(ecat63errmsg, "matrix list is empty"); fclose(fp); return(5);}
614  /* Sort matrix list */
615  for(i=0; i<mlist.matrixNr-1; i++) for(j=i+1; j<mlist.matrixNr; j++) {
616  if(mlist.matdir[i].matnum>mlist.matdir[j].matnum) {
617  tmpMatdir=mlist.matdir[i];
618  mlist.matdir[i]=mlist.matdir[j]; mlist.matdir[j]=tmpMatdir;
619  }
620  }
621  /* Trim extra frames */
622  if(main_header.num_frames>0) {
623  del_nr=ecat63DeleteLateFrames(&mlist, main_header.num_frames);
624  if(IMG_TEST && del_nr>0)
625  printf(" %d entries in matrix list will not be used.\n", del_nr);
626  }
627 
628  /* Decide the plane to read */
629  if(img->status==IMG_STATUS_OCCUPIED) {
630  /* img contains data */
631  isFirst=0; prev_frameNr=img->dimt;
632  /* get the current plane in there */
633  prev_plane=img->planeNumber[img->dimz-1];
634  /* Clear all data in img: not here but only after finding new data */
635  } else {
636  /* img does not contain any data */
637  isFirst=1;
638  /* set "current plane" to -1 */
639  prev_plane=-1;
640  }
641  /* from sorted matrix list, find the first plane larger than the current one */
642  for(m=0, plane=-1; m<mlist.matrixNr; m++) if(mlist.matdir[m].matstat==1) {
643  mat_numdoc(mlist.matdir[m].matnum, &matval);
644  if(matval.plane>prev_plane) {plane=matval.plane; break;}
645  }
646  /* If not found, return an error or 1 if this was not the first one */
647  if(plane<0) {
648  fclose(fp); ecat63EmptyMatlist(&mlist);
649  if(isFirst) {
650  sprintf(ecat63errmsg, "ECAT file contains no matrices");
651  return(7);
652  } else {
653  sprintf(ecat63errmsg, "ECAT file contains no more planes");
654  if(IMG_TEST) printf("%s\n", ecat63errmsg);
655  return(1);
656  }
657  }
658  if(IMG_TEST) printf("Next plane to read: %d\n", plane);
659  /* Clear all data in img */
660  imgEmpty(img);
661 
662  /* Calculate the number of frames and gates */
663  prev_frame=frame=-1; frameNr=0; ret=0; planeNr=1;
664  for(m=0; m<mlist.matrixNr; m++) if(mlist.matdir[m].matstat==1) {
665  /* get frame and plane; work only with current plane */
666  mat_numdoc(mlist.matdir[m].matnum, &matval);
667  if(matval.plane<plane) continue; else if(matval.plane>plane) break;
668  if(main_header.num_frames>=main_header.num_gates) frame=matval.frame;
669  else frame=matval.gate;
670  frameNr++;
671  /* frame (gate) numbers must be continuous */
672  if(prev_frame>0 && frame!=prev_frame+1) {ret=1; break;}
673  prev_frame=frame;
674  /* Calculate and check the size of one data matrix */
675  if(frameNr==1) {
676  firstm=m;
677  blkNr=mlist.matdir[m].endblk-mlist.matdir[m].strtblk;
678  } else if(blkNr!=mlist.matdir[m].endblk-mlist.matdir[m].strtblk) {
679  ret=2; break;
680  }
681  prev_frame=frame;
682  } /* next matrix */
683  if(ret==1) {
684  strcpy(ecat63errmsg, "missing matrix");
685  ecat63EmptyMatlist(&mlist); fclose(fp); return(6);
686  }
687  if(ret==2) {
688  strcpy(ecat63errmsg, "matrix sizes are different");
689  ecat63EmptyMatlist(&mlist); fclose(fp); return(6);
690  }
691  /* Check that frameNr is equal to the dimt in IMG */
692  if(!isFirst && frameNr!=prev_frameNr) {
693  strcpy(ecat63errmsg, "different frame nr in different planes");
694  ecat63EmptyMatlist(&mlist); fclose(fp); return(6);
695  }
696 
697  /* Read x,y-dimensions from 1st matrix subheader on current plane */
698  m=firstm; if(main_header.file_type==IMAGE_DATA) {
699  ret=ecat63ReadImageheader(fp, mlist.matdir[m].strtblk, &image_header);
700  dim_x=image_header.dimension_1; dim_y=image_header.dimension_2;
701  } else if(main_header.file_type==RAW_DATA) {
702  ret=ecat63ReadScanheader(fp, mlist.matdir[m].strtblk, &scan_header);
703  dim_x=scan_header.dimension_1; dim_y=scan_header.dimension_2;
704  } else if(main_header.file_type==ATTN_DATA) {
705  ret=ecat63ReadAttnheader(fp, mlist.matdir[m].strtblk, &attn_header);
706  dim_x=attn_header.dimension_1; dim_y=attn_header.dimension_2;
707  } else if(main_header.file_type==NORM_DATA) {
708  ret=ecat63ReadNormheader(fp, mlist.matdir[m].strtblk, &norm_header);
709  dim_x=norm_header.dimension_1; dim_y=norm_header.dimension_2;
710  }
711  pxlNr=dim_x*dim_y;
712  if(ret) {
713  sprintf(ecat63errmsg, "cannot read matrix %u subheader", mlist.matdir[m].matnum);
714  ecat63EmptyMatlist(&mlist); fclose(fp); return(7);
715  }
716 
717  /* Allocate memory for ECAT data matrix */
718  if(IMG_TEST) printf("allocating memory for %d blocks\n", blkNr);
719  mdata=(char*)malloc(blkNr*MatBLKSIZE); if(mdata==NULL) {
720  strcpy(ecat63errmsg, "out of memory");
721  fclose(fp); ecat63EmptyMatlist(&mlist); return(8);
722  }
723  /* Allocate memory for whole img data */
724  ret=imgAllocate(img, planeNr, dim_y, dim_x, frameNr);
725  if(ret) {
726  sprintf(ecat63errmsg, "out of memory (%d)", ret);
727  fclose(fp); ecat63EmptyMatlist(&mlist); free(mdata); return(8);
728  }
729 
730  /* Fill img info with ECAT main and sub header information */
731  img->scanner=main_header.system_type;
732  img->unit=main_header.calibration_units; /* this continues below */
733  strncpy(img->radiopharmaceutical, main_header.radiopharmaceutical, 32);
734  img->isotopeHalflife=main_header.isotope_halflife;
735  memset(&scanStart, 0, sizeof(struct tm));
736  scanStart.tm_year=main_header.scan_start_year-1900;
737  scanStart.tm_mon=main_header.scan_start_month-1;
738  scanStart.tm_mday=main_header.scan_start_day; scanStart.tm_yday=0;
739  scanStart.tm_hour=main_header.scan_start_hour;
740  scanStart.tm_min=main_header.scan_start_minute;
741  scanStart.tm_sec=main_header.scan_start_second;
742  scanStart.tm_isdst=-1;
743  img->scanStart=mktime(&scanStart); /*printf("img->scanStart=%d\n", img->scanStart);*/
744  if(img->scanStart==-1) img->scanStart=0;
745  img->axialFOV=10.*main_header.axial_fov;
746  img->transaxialFOV=10.*main_header.transaxial_fov;
747  strncpy(img->studyNr, main_header.study_name, MAX_STUDYNR_LEN);
748  strncpy(img->patientName, main_header.patient_name, 31);
749  strncpy(img->patientID, main_header.patient_id, 15);
750  img->sizez=10.*main_header.plane_separation;
751  if(main_header.file_type==IMAGE_DATA) {
752  img->type=IMG_TYPE_IMAGE;
753  if(img->unit<1) img->unit=image_header.quant_units;
754  img->zoom=image_header.recon_scale;
755  if(image_header.decay_corr_fctr>1.0) img->decayCorrected=1;
756  img->sizex=img->sizey=10.*image_header.pixel_size;
757  if(img->sizez<10.*image_header.slice_width)
758  img->sizez=10.*image_header.slice_width;
759  } else if(main_header.file_type==RAW_DATA) {
760  img->type=IMG_TYPE_RAW;
761  img->decayCorrected=0;
762  } else if(main_header.file_type==ATTN_DATA) {
763  img->type=IMG_TYPE_RAW;
764  img->decayCorrected=0;
765  } else if(main_header.file_type==NORM_DATA) {
766  img->type=IMG_TYPE_RAW;
767  img->decayCorrected=0;
768  }
769  img->planeNumber[0]=plane;
770  strncpy(img->studyDescription, main_header.study_description, 32);
771  strncpy(img->userProcessCode, main_header.user_process_code, 10); img->userProcessCode[10]=(char)0;
772  /* If valid study number is found in user_process_code, then take it */
773  if(!img->studyNr[0] && studynr_validity_check(img->userProcessCode))
774  strcpy(img->studyNr, img->userProcessCode);
775  /* Set _fileFormat */
776  img->_fileFormat=IMG_E63;
777 
778  /* Read one ECAT matrix at a time and put them to img */
779  for(m=firstm, frame=1; m<mlist.matrixNr; m++) if(mlist.matdir[m].matstat==1) {
780  /* get plane */
781  mat_numdoc(mlist.matdir[m].matnum, &matval);
782  if(matval.plane>plane) break; /* Quit when current plane is read */
783  /* Read subheader */
784  if(main_header.file_type==IMAGE_DATA) {
785  ret=ecat63ReadImageheader(fp, mlist.matdir[m].strtblk, &image_header);
786  if(dim_x!=image_header.dimension_1 || dim_y!=image_header.dimension_2) ret=1;
787  scale=image_header.quant_scale;
788  if(image_header.ecat_calibration_fctr>0.0
789  && fabs(main_header.calibration_factor/image_header.ecat_calibration_fctr-1.0)>0.0001)
790  scale*=image_header.ecat_calibration_fctr;
791  if(img->unit==0 && image_header.quant_units>0) img->unit=image_header.quant_units;
792  datatype=image_header.data_type;
793  img->start[frame-1]=image_header.frame_start_time/1000.;
794  img->end[frame-1]=img->start[frame-1]+image_header.frame_duration/1000.;
795  img->mid[frame-1]=0.5*(img->start[frame-1]+img->end[frame-1]);
796  if(image_header.decay_corr_fctr>1.0)
797  img->decayCorrFactor[frame-1]=image_header.decay_corr_fctr;
798  } else if(main_header.file_type==RAW_DATA) {
799  ret=ecat63ReadScanheader(fp, mlist.matdir[m].strtblk, &scan_header);
800  if(dim_x!=scan_header.dimension_1 || dim_y!=scan_header.dimension_2) ret=1;
801  scale=scan_header.scale_factor;
802  if(scan_header.loss_correction_fctr>0.0) scale*=scan_header.loss_correction_fctr;
803  datatype=scan_header.data_type;
804  img->start[frame-1]=scan_header.frame_start_time/1000.;
805  img->end[frame-1]=img->start[frame-1]+scan_header.frame_duration/1000.;
806  img->mid[frame-1]=0.5*(img->start[frame-1]+img->end[frame-1]);
807  img->sampleDistance=10.0*scan_header.sample_distance;
808  img->prompts[frame-1]=(float)scan_header.prompts;
809  img->randoms[frame-1]=(float)scan_header.delayed;
810  } else if(main_header.file_type==ATTN_DATA) {
811  ret=ecat63ReadAttnheader(fp, mlist.matdir[m].strtblk, &attn_header);
812  if(dim_x!=attn_header.dimension_1 || dim_y!=attn_header.dimension_2) ret=1;
813  scale=attn_header.scale_factor;
814  datatype=attn_header.data_type;
815  } else if(main_header.file_type==NORM_DATA) {
816  ret=ecat63ReadNormheader(fp, mlist.matdir[m].strtblk, &norm_header);
817  if(dim_x!=norm_header.dimension_1 || dim_y!=norm_header.dimension_2) ret=1;
818  scale=norm_header.scale_factor;
819  datatype=norm_header.data_type;
820  } else datatype=-1;
821  if(ret) {
822  sprintf(ecat63errmsg, "cannot read matrix %u subheader", mlist.matdir[m].matnum);
823  ecat63EmptyMatlist(&mlist); free(mdata); fclose(fp); return(7);
824  }
825  if(main_header.calibration_factor>0.0) scale*=main_header.calibration_factor;
826  if(IMG_TEST>2) printf("scale=%e datatype=%d\n", scale, datatype);
827  /* Read the pixel data */
828  ret=ecat63ReadMatdata(fp, mlist.matdir[m].strtblk+1,
829  mlist.matdir[m].endblk-mlist.matdir[m].strtblk,
830  mdata, datatype);
831  if(ret) {
832  strcpy(ecat63errmsg, "cannot read matrix data");
833  ecat63EmptyMatlist(&mlist); free(mdata); fclose(fp); return(9);
834  }
835  if(datatype==BYTE_TYPE) {
836  for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr++) {
837  img->m[0][i][j][frame-1]=scale*(float)(*mptr);
838  }
839  } else if(datatype==VAX_I2 || datatype==SUN_I2) {
840  for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr+=2) {
841  sptr=(short int*)mptr;
842  img->m[0][i][j][frame-1]=scale*(float)(*sptr);
843  }
844  } else if(datatype==VAX_I4 || datatype==SUN_I4) {
845  for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr+=4) {
846  iptr=(int*)mptr;
847  img->m[0][i][j][frame-1]=scale*(float)(*iptr);
848  }
849  } else if(datatype==VAX_R4 || datatype==IEEE_R4) {
850  for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr+=4) {
851  memcpy(&img->m[0][i][j][frame-1], mptr, 4);
852  img->m[0][i][j][frame-1]*=scale;
853  }
854  }
855  frame++;
856  } /* next matrix (frame) */
857  /* Set the unit */
858  imgUnitFromEcat(img, img->unit);
859 
860  /* Free memory and close file */
861  free(mdata);
862  ecat63EmptyMatlist(&mlist);
863  fclose(fp);
864 
865  /* Set _dataType if it has not been set before */
866  if(img->_dataType<1) img->_dataType=datatype;
867  /* For saving, only 2-byte VAX data types are allowed */
868  if(img->_dataType==VAX_I4 || img->_dataType==VAX_R4) img->_dataType=VAX_I2;
869 
870  return(0);
871 }
872 /*****************************************************************************/
873 
874 /*****************************************************************************/
886 int ecat63AddImg(const char *fname, IMG *img) {
887  int n, i, j, m, ret=0, add=0;
888  int frameNr, planeNr;
889  int frame, plane, prev_frame, prev_plane;
890  float f, fmax, fmin, g, scale;
891  short int *sdata, *sptr, smin, smax;
892  FILE *fp;
893  ECAT63_mainheader main_header;
894  ECAT63_imageheader image_header;
895  ECAT63_scanheader scan_header;
896  MATRIXLIST mlist;
897  MatDir tmpMatdir;
898  Matval matval;
899  struct tm *scanStart;
900 
901 
902  if(IMG_TEST) printf("ecat63AddImg(%s, *img)\n", fname);
903  if(IMG_TEST>1) imgInfo(img);
904  /* Check the arguments */
905  if(fname==NULL) {strcpy(ecat63errmsg, "invalid ECAT filename"); return(1);}
906  if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) {
907  strcpy(ecat63errmsg, "image data is empty"); return(2);}
908  if(img->_dataType<1) img->_dataType=VAX_I2;
909 
910  /* Initiate headers */
911  memset(&main_header, 0, sizeof(ECAT63_mainheader));
912  memset(&image_header,0, sizeof(ECAT63_imageheader));
913  memset(&scan_header, 0, sizeof(ECAT63_scanheader));
914 
915  /* Make a main header */
916  main_header.sw_version=2;
917  main_header.num_planes=img->dimz;
918  main_header.num_frames=img->dimt;
919  main_header.num_gates=1;
920  main_header.num_bed_pos=1;
921  if(img->type==IMG_TYPE_IMAGE) main_header.file_type=IMAGE_DATA;
922  else main_header.file_type=RAW_DATA;
923  main_header.data_type=img->_dataType;
924  if(img->scanner>0) main_header.system_type=img->scanner;
925  else main_header.system_type=ECAT63_SYSTEM_TYPE_DEFAULT;
926  main_header.calibration_factor=1.0;
927  main_header.axial_fov=img->axialFOV/10.0;
928  main_header.transaxial_fov=img->transaxialFOV/10.0;
929  main_header.plane_separation=img->sizez/10.0;
930  main_header.calibration_units=imgUnitToEcat6(img);
931  strncpy(main_header.radiopharmaceutical, img->radiopharmaceutical, 32);
932  scanStart=localtime(&img->scanStart);
933  if(scanStart!=NULL) {
934  main_header.scan_start_year=scanStart->tm_year+1900;
935  main_header.scan_start_month=scanStart->tm_mon+1;
936  main_header.scan_start_day=scanStart->tm_mday;
937  main_header.scan_start_hour=scanStart->tm_hour;
938  main_header.scan_start_minute=scanStart->tm_min;
939  main_header.scan_start_second=scanStart->tm_sec;
940  } else {
941  main_header.scan_start_year=1900;
942  main_header.scan_start_month=1;
943  main_header.scan_start_day=1;
944  main_header.scan_start_hour=0;
945  main_header.scan_start_minute=0;
946  main_header.scan_start_second=0;
947  }
948  main_header.isotope_halflife=img->isotopeHalflife;
949  strcpy(main_header.isotope_code, imgIsotope(img));
950  strcpy(main_header.study_name, img->studyNr);
951  strcpy(main_header.patient_name, img->patientName);
952  strcpy(main_header.patient_id, img->patientID);
953  strncpy(main_header.study_description, img->studyDescription, 32);
954  strncpy(main_header.user_process_code, img->userProcessCode, 10);
955 
956  /* Check if the ECAT file exists */
957  if(access(fname, 0) != -1) {
958  if(IMG_TEST) printf("Opening existing ECAT file.\n");
959  add=1;
960  /* Open the ECAT file */
961  if((fp=fopen(fname, "rb+")) == NULL) {
962  strcpy(ecat63errmsg, "cannot create ECAT file"); return(3);}
963  /* Read main header */
964  if((ret=ecat63ReadMainheader(fp, &main_header))) {
965  sprintf(ecat63errmsg, "cannot read main header (%d)", ret);
966  fclose(fp); return(3);
967  }
968  fflush(fp);
969  /* Check that filetypes are matching */
970  if((main_header.file_type==IMAGE_DATA && img->type==IMG_TYPE_IMAGE) ||
971  (main_header.file_type==RAW_DATA && img->type==IMG_TYPE_RAW)) {
972  /* ok */
973  } else {
974  sprintf(ecat63errmsg, "cannot add different filetype");
975  fclose(fp); return(3);
976  }
977  } else {
978  if(IMG_TEST) printf("ECAT file does not exist.\n");
979  add=0;
980  /* Create output ECAT file */
981  fp=ecat63Create(fname, &main_header);
982  if(fp==NULL) {strcpy(ecat63errmsg, "cannot create ECAT file"); return(3);}
983  }
984  if(IMG_TEST) ecat63PrintMainheader(&main_header, stdout);
985 
986  /* Allocate memory for matrix data */
987  sdata=(short int*)malloc(img->dimx*img->dimy*sizeof(short int));
988  if(sdata==NULL) {strcpy(ecat63errmsg, "out of memory"); return(4);}
989 
990  /* Set the subheader contents, as far as possible */
991  switch(main_header.file_type) {
992  case RAW_DATA:
993  scan_header.data_type=main_header.data_type;
994  scan_header.dimension_1=img->dimx;
995  scan_header.dimension_2=img->dimy;
996  scan_header.frame_duration_sec=1;
997  scan_header.scale_factor=1.0;
998  scan_header.frame_start_time=0;
999  scan_header.frame_duration=1000;
1000  scan_header.loss_correction_fctr=1.0;
1001  scan_header.sample_distance=(img->sampleDistance)/10.0;
1002  /*if(IMG_TEST) ecat63PrintScanheader(&scan_header);*/
1003  break;
1004  case IMAGE_DATA:
1005  image_header.data_type=main_header.data_type;
1006  image_header.num_dimensions=2;
1007  image_header.dimension_1=img->dimx;
1008  image_header.dimension_2=img->dimy;
1009  image_header.recon_scale=img->zoom;
1010  image_header.quant_scale=1.0;
1011  image_header.slice_width=img->sizez/10.;
1012  image_header.pixel_size=img->sizex/10.;
1013  image_header.frame_start_time=0;
1014  image_header.frame_duration=1000;
1015  image_header.plane_eff_corr_fctr=1.0;
1016  image_header.decay_corr_fctr=0.0;
1017  image_header.loss_corr_fctr=1.0;
1018  image_header.ecat_calibration_fctr=1.0;
1019  image_header.well_counter_cal_fctr=1.0;
1020  image_header.quant_units=main_header.calibration_units;
1021  /*if(IMG_TEST) ecat63PrintImageheader(&image_header);*/
1022  break;
1023  }
1024 
1025  /* Write one matrix at a time */
1026  n=0;
1027  for(plane=1; plane<=img->dimz; plane++) for(frame=1; frame<=img->dimt; frame++) {
1028  /* Scale data to short ints */
1029  /* Search min and max */
1030  fmin=fmax=f=img->m[plane-1][0][0][frame-1];
1031  for(i=0; i<img->dimy; i++) for(j=0; j<img->dimx; j++) {
1032  f=img->m[plane-1][i][j][frame-1];
1033  if(f>fmax) fmax=f; else if(f<fmin) fmin=f;
1034  }
1035  /* Calculate scaling factor */
1036  if(fabs(fmin)>fabs(fmax)) g=fabs(fmin); else g=fabs(fmax);
1037  if(g!=0) scale=32766./g; else scale=1.0;
1038  /*printf("fmin=%e fmax=%e g=%e scale=%e\n", fmin, fmax, g, scale);*/
1039  /* Scale matrix data to shorts */
1040  sptr=sdata;
1041  for(i=0; i<img->dimy; i++) for(j=0; j<img->dimx; j++) {
1042  *sptr=(short int)temp_roundf(scale*img->m[plane-1][i][j][frame-1]);
1043  sptr++;
1044  }
1045  /* Calculate and set subheader min&max and scale */
1046  smin=(short int)temp_roundf(scale*fmin); smax=(short int)temp_roundf(scale*fmax);
1047  if(main_header.file_type==RAW_DATA) {
1048  scan_header.scan_min=smin; scan_header.scan_max=smax;
1049  scan_header.scale_factor=1.0/scale;
1050  scan_header.frame_start_time=(int)temp_roundf(1000.*img->start[frame-1]);
1051  scan_header.frame_duration=
1052  (int)temp_roundf(1000.*(img->end[frame-1]-img->start[frame-1]));
1053  scan_header.prompts=temp_roundf(img->prompts[frame-1]);
1054  scan_header.delayed=temp_roundf(img->randoms[frame-1]);
1055  } else if(main_header.file_type==IMAGE_DATA) {
1056  image_header.image_min=smin; image_header.image_max=smax;
1057  image_header.quant_scale=1.0/scale;
1058  image_header.frame_start_time=(int)temp_roundf(1000.*img->start[frame-1]);
1059  image_header.frame_duration=
1060  (int)temp_roundf(1000.*(img->end[frame-1]-img->start[frame-1]));
1061  /* Set decay correction factor */
1062  if(img->decayCorrected)
1063  image_header.decay_corr_fctr=img->decayCorrFactor[frame-1];
1064  }
1065  /* Write matrix data */
1066  m=mat_numcod(frame, img->planeNumber[plane-1], 1, 0, 0);
1067  sptr=sdata;
1068  if(IMG_TEST) printf(" writing matnum=%d\n", m);
1069  if(main_header.file_type==RAW_DATA) {
1070  /*if(IMG_TEST) ecat63PrintScanheader(&scan_header);*/
1071  ret=ecat63WriteScan(fp, m, &scan_header, sptr);
1072  } else if(main_header.file_type==IMAGE_DATA) {
1073  /*if(IMG_TEST) ecat63PrintImageheader(&image_header);*/
1074  ret=ecat63WriteImage(fp, m, &image_header, sptr);
1075  }
1076  if(ret) {
1077  sprintf(ecat63errmsg, "cannot write data on pl%02d fr%02d (%d).",
1078  plane, frame, ret);
1079  fclose(fp); remove(fname); free(sdata);
1080  return(9);
1081  }
1082  n++;
1083  } /* next matrix */
1084  if(IMG_TEST) printf(" %d matrices written in %s\n", n, fname);
1085 
1086  /* Free matrix memory */
1087  free(sdata);
1088 
1089  /* If matrices were added, set main header */
1090  if(add==1) {
1091  fflush(fp);
1092  /* Read matrix list */
1093  ecat63InitMatlist(&mlist);
1094  ret=ecat63ReadMatlist(fp, &mlist);
1095  if(ret) {
1096  sprintf(ecat63errmsg, "cannot read matrix list (%d)", ret);
1097  fclose(fp); return(21);
1098  }
1099  if(mlist.matrixNr<=0) {
1100  strcpy(ecat63errmsg, "matrix list is empty"); fclose(fp); return(21);}
1101  /* Sort matrix list */
1102  for(i=0; i<mlist.matrixNr-1; i++) for(j=i+1; j<mlist.matrixNr; j++) {
1103  if(mlist.matdir[i].matnum>mlist.matdir[j].matnum) {
1104  tmpMatdir=mlist.matdir[i];
1105  mlist.matdir[i]=mlist.matdir[j]; mlist.matdir[j]=tmpMatdir;
1106  }
1107  }
1108  /* Calculate the number of planes and frames */
1109  prev_plane=plane=-1; prev_frame=frame=-1; frameNr=planeNr=0;
1110  for(m=0; m<mlist.matrixNr; m++) {
1111  mat_numdoc(mlist.matdir[m].matnum, &matval);
1112  plane=matval.plane; frame=matval.frame;
1113  if(plane!=prev_plane) {frameNr=1; planeNr++;} else {frameNr++;}
1114  prev_plane=plane; prev_frame=frame;
1115  } /* next matrix */
1116  ecat63EmptyMatlist(&mlist);
1117  main_header.num_planes=planeNr; main_header.num_frames=frameNr;
1118  /* Write main header */
1119  ret=ecat63WriteMainheader(fp, &main_header);
1120  if(ret) {
1121  sprintf(ecat63errmsg, "cannot write mainheader (%d)", ret);
1122  fclose(fp); return(22);
1123  }
1124  }
1125 
1126  /* Close file and free memory */
1127  fclose(fp);
1128 
1129  return(0);
1130 }
1131 /*****************************************************************************/
1132 
1133 /*****************************************************************************/
1141  if(h==NULL) return(0);
1142  if(h->file_type==IMAGE_DATA) return(1);
1143  if(h->file_type==RAW_DATA) return(1);
1144  if(h->file_type==ATTN_DATA) return(1);
1145  if(h->file_type==NORM_DATA) return(1);
1146  return(0);
1147 }
1148 /*****************************************************************************/
1149 
1150 /*****************************************************************************/
1158  struct tm scanStart;
1159 
1160  img->_dataType=h->data_type; /* again in subheaders*/
1161  /* For saving IMG data, only 2-byte VAX data types are allowed, so change it now */
1162  if(img->_dataType==VAX_I4 || img->_dataType==VAX_R4) img->_dataType=VAX_I2;
1163  img->scanner=h->system_type;
1164  strncpy(img->radiopharmaceutical, h->radiopharmaceutical, 32);
1166  {
1167  memset(&scanStart, 0, sizeof(struct tm));
1168  scanStart.tm_year=h->scan_start_year-1900;
1169  scanStart.tm_mon=h->scan_start_month-1;
1170  scanStart.tm_mday=h->scan_start_day; scanStart.tm_yday=0;
1171  scanStart.tm_hour=h->scan_start_hour;
1172  scanStart.tm_min=h->scan_start_minute;
1173  scanStart.tm_sec=h->scan_start_second;
1174  scanStart.tm_isdst=-1;
1175  img->scanStart=mktime(&scanStart); if(img->scanStart==-1) img->scanStart=0;
1176  }
1177  img->axialFOV=10.0*h->axial_fov;
1178  img->transaxialFOV=10.0*h->transaxial_fov;
1179  strncpy(img->studyNr, h->study_name, MAX_STUDYNR_LEN);
1180  strncpy(img->patientName, h->patient_name, 31);
1181  strncpy(img->patientID, h->patient_id, 15);
1182  img->sizez=10.0*h->plane_separation;
1183  switch(h->file_type) {
1184  case IMAGE_DATA: img->type=IMG_TYPE_IMAGE; break;
1185  case RAW_DATA:
1186  case ATTN_DATA:
1187  case NORM_DATA:
1188  default: img->type=IMG_TYPE_RAW;
1189  }
1190  strncpy(img->studyDescription, h->study_description, 32);
1191  strncpy(img->userProcessCode, h->user_process_code, 10);
1192  img->userProcessCode[10]=(char)0;
1193  /* If valid study number is found in user_process_code, then take it */
1194  if(!img->studyNr[0] && studynr_validity_check(img->userProcessCode))
1195  strcpy(img->studyNr, img->userProcessCode);
1196  /* Set calibration unit (this may have to be read from subheader later) */
1198 }
1199 /*****************************************************************************/
1200 
1201 /*****************************************************************************/
1209  struct tm *scanStart;
1210 
1211  h->sw_version=2;
1212  h->num_planes=img->dimz;
1213  h->num_frames=img->dimt;
1214  h->num_gates=1;
1215  h->num_bed_pos=1;
1216  if(img->type==IMG_TYPE_IMAGE) h->file_type=IMAGE_DATA;
1217  else h->file_type=RAW_DATA;
1218  h->data_type=VAX_I2;
1219  if(img->scanner>0) h->system_type=img->scanner;
1221  h->calibration_factor=1.0;
1222  h->axial_fov=img->axialFOV/10.0;
1223  h->transaxial_fov=img->transaxialFOV/10.0;
1224  h->plane_separation=img->sizez/10.0;
1226  strncpy(h->radiopharmaceutical, img->radiopharmaceutical, 32);
1227  scanStart=localtime(&img->scanStart);
1228  if(scanStart!=NULL) {
1229  h->scan_start_year=scanStart->tm_year+1900;
1230  h->scan_start_month=scanStart->tm_mon+1;
1231  h->scan_start_day=scanStart->tm_mday;
1232  h->scan_start_hour=scanStart->tm_hour;
1233  h->scan_start_minute=scanStart->tm_min;
1234  h->scan_start_second=scanStart->tm_sec;
1235  } else {
1236  h->scan_start_year=1900;
1237  h->scan_start_month=1;
1238  h->scan_start_day=1;
1239  h->scan_start_hour=0;
1240  h->scan_start_minute=0;
1241  h->scan_start_second=0;
1242  }
1244  strcpy(h->isotope_code, imgIsotope(img));
1245  strcpy(h->study_name, img->studyNr);
1246  strcpy(h->patient_name, img->patientName);
1247  strcpy(h->patient_id, img->patientID);
1248  strncpy(h->user_process_code, img->userProcessCode, 10);
1249  strncpy(h->study_description, img->studyDescription, 32);
1250 }
1251 /*****************************************************************************/
1252 
1253 /*****************************************************************************/
1261  int fileFormat=IMG_UNKNOWN;
1262  switch(h->file_type) {
1263  case IMAGE_DATA:
1264  case RAW_DATA:
1265  case ATTN_DATA:
1266  case NORM_DATA:
1267  fileFormat=IMG_E63; break;
1268  default:
1269  fileFormat=IMG_UNKNOWN; break;
1270  }
1271  return fileFormat;
1272 }
1273 /*****************************************************************************/
1274 
1275 /*****************************************************************************/
1289 int imgReadEcat63Header(const char *fname, IMG *img) {
1290  int m, blkNr=0, ret=0;
1291  int frameNr, planeNr, del_nr=0;
1292  FILE *fp;
1293  ECAT63_mainheader main_header;
1294  MATRIXLIST mlist;
1295  ECAT63_imageheader image_header;
1296  ECAT63_scanheader scan_header;
1297  ECAT63_attnheader attn_header;
1298  ECAT63_normheader norm_header;
1299 
1300 
1301  if(IMG_TEST) printf("\nimgReadEcat63Header(%s, *img)\n", fname);
1302 
1303  /* Check the arguments */
1304  if(img==NULL) return STATUS_FAULT;
1305  if(img->status!=IMG_STATUS_INITIALIZED) return STATUS_FAULT;
1306  imgSetStatus(img, STATUS_FAULT);
1307  if(fname==NULL) return STATUS_FAULT;
1308 
1309  /* Open the file */
1310  if((fp=fopen(fname, "rb")) == NULL) return STATUS_NOFILE;
1311 
1312  /* Read main header */
1313  if((ret=ecat63ReadMainheader(fp, &main_header))) {
1314  fclose(fp); return STATUS_NOMAINHEADER;}
1315  /* Check if file_type is supported */
1316  if(imgEcat63Supported(&main_header)==0) {fclose(fp); return STATUS_UNSUPPORTED;}
1317  /* Copy main header information into IMG; sets also img.type */
1318  imgGetEcat63MHeader(img, &main_header);
1319  if(IMG_TEST>7) printf("img.type := %d\n", img->type);
1320  img->_fileFormat=imgGetEcat63Fileformat(&main_header);
1321  if(IMG_TEST>7) printf("img._fileFormat := %d\n", img->_fileFormat);
1322  if(img->_fileFormat==IMG_UNKNOWN) {fclose(fp); return STATUS_UNSUPPORTED;}
1323 
1324  /* Read matrix list and nr and sort it */
1325  ecat63InitMatlist(&mlist);
1326  ret=ecat63ReadMatlist(fp, &mlist);
1327  if(ret) {fclose(fp); return STATUS_NOMATLIST;}
1328  if(mlist.matrixNr<=0) {
1329  ecat63EmptyMatlist(&mlist); fclose(fp); return STATUS_INVALIDMATLIST;}
1330  /* Make sure that plane and frame numbers are continuous */
1331  ecat63GatherMatlist(&mlist, 1, 1, 1, 1);
1332  ecat63SortMatlistByPlane(&mlist); /* otherwise frameNr cannot be computed as below */
1333  /* Trim extra frames */
1334  if(main_header.num_frames>0)
1335  del_nr=ecat63DeleteLateFrames(&mlist, main_header.num_frames);
1336  /* Get plane and frame numbers and check that volume is full */
1337  ret=ecat63GetPlaneAndFrameNr(&mlist, &main_header, &planeNr, &frameNr);
1338  if(ret) {ecat63EmptyMatlist(&mlist); fclose(fp); return ret;}
1339  img->dimz=planeNr;
1340  img->dimt=frameNr;
1341  /* Calculate and check the size of one data matrix */
1342  ret=ecat63GetMatrixBlockSize(&mlist, &blkNr);
1343  if(ret) {ecat63EmptyMatlist(&mlist); fclose(fp); return ret;}
1344 
1345  /* Read one subheader */
1346  if(IMG_TEST>5) printf("main_header.file_type := %d\n", main_header.file_type);
1347  m=0;
1348  switch(main_header.file_type) {
1349  case IMAGE_DATA:
1350  ret=ecat63ReadImageheader(fp, mlist.matdir[m].strtblk, &image_header);
1351  break;
1352  case RAW_DATA:
1353  ret=ecat63ReadScanheader(fp, mlist.matdir[m].strtblk, &scan_header);
1354  break;
1355  case ATTN_DATA:
1356  ret=ecat63ReadAttnheader(fp, mlist.matdir[m].strtblk, &attn_header);
1357  break;
1358  case NORM_DATA:
1359  ret=ecat63ReadNormheader(fp, mlist.matdir[m].strtblk, &norm_header);
1360  break;
1361  default: ret=-1;
1362  }
1363  /* Free locally allocated memory and close the file */
1364  ecat63EmptyMatlist(&mlist); fclose(fp);
1365  /* Check whether subheader was read */
1366  if(ret) return STATUS_NOSUBHEADER;
1367 
1368  /* Get the following information in the subheader:
1369  dimensions x, y and z; datatype;
1370  image decay correction on/off, zoom, and pixel size;
1371  sinogram sample distance.
1372  */
1373  switch(main_header.file_type) {
1374  case IMAGE_DATA:
1375  img->dimx=image_header.dimension_1; img->dimy=image_header.dimension_2;
1376  if(img->unit<1) imgUnitFromEcat(img, image_header.quant_units);
1377  img->_dataType=image_header.data_type;
1378  img->zoom=image_header.recon_scale;
1379  if(image_header.decay_corr_fctr>1.0) img->decayCorrected=1;
1380  img->sizex=img->sizey=10.*image_header.pixel_size;
1381  if(img->sizez<10.*image_header.slice_width)
1382  img->sizez=10.*image_header.slice_width;
1383  break;
1384  case RAW_DATA:
1385  img->dimx=scan_header.dimension_1; img->dimy=scan_header.dimension_2;
1386  img->type=IMG_TYPE_RAW;
1387  img->_dataType=scan_header.data_type;
1388  img->decayCorrected=0;
1389  img->sampleDistance=10.0*scan_header.sample_distance;
1390  break;
1391  case ATTN_DATA:
1392  img->dimx=attn_header.dimension_1; img->dimy=attn_header.dimension_2;
1393  img->type=IMG_TYPE_RAW;
1394  img->decayCorrected=0;
1395  img->_dataType=attn_header.data_type;
1396  break;
1397  case NORM_DATA:
1398  img->dimx=norm_header.dimension_1; img->dimy=norm_header.dimension_2;
1399  img->type=IMG_TYPE_RAW;
1400  img->decayCorrected=0;
1401  img->_dataType=norm_header.data_type;
1402  break;
1403  }
1404 
1405  /* For saving IMG data, only 2-byte VAX data types are allowed, so change it now */
1406  if(img->_dataType==VAX_I4 || img->_dataType==VAX_R4) img->_dataType=VAX_I2;
1407 
1408  imgSetStatus(img, STATUS_OK);
1409  return STATUS_OK;
1410 }
1411 /*****************************************************************************/
1412 
1413 /*****************************************************************************/
1422 int imgReadEcat63FirstFrame(const char *fname, IMG *img) {
1423  int ret=0;
1424 
1425  if(IMG_TEST) printf("\nimgReadEcat63FirstFrame(%s, *img)\n", fname);
1426  /* Check the input */
1427  if(img==NULL) return STATUS_FAULT;
1428  if(img->status!=IMG_STATUS_INITIALIZED) return STATUS_FAULT;
1429  imgSetStatus(img, STATUS_FAULT);
1430  if(fname==NULL) return STATUS_FAULT;
1431 
1432  /* Read header information from file */
1433  ret=imgReadEcat63Header(fname, img); if(ret) return(ret);
1434  if(IMG_TEST>3) imgInfo(img);
1435 
1436  /* Allocate memory for one frame */
1437  img->dimt=1;
1438  ret=imgAllocate(img, img->dimz, img->dimy, img->dimx, img->dimt);
1439  if(ret) return STATUS_NOMEMORY;
1440 
1441  /* Read the first frame */
1442  ret=imgReadEcat63Frame(fname, 1, img, 0); if(ret) return(ret);
1443 
1444  /* All went well */
1445  imgSetStatus(img, STATUS_OK);
1446  return STATUS_OK;
1447 }
1448 /*****************************************************************************/
1449 
1450 /*****************************************************************************/
1464 int imgReadEcat63Frame(const char *fname, int frame_to_read, IMG *img, int frame_index) {
1465  int m, ret=0, blkNr=0, frame, plane, seqplane=-1, xi, yi;
1466  int local_data_type=0;
1467  FILE *fp;
1468  ECAT63_mainheader main_header;
1469  MATRIXLIST mlist;
1470  Matval matval;
1471  ECAT63_imageheader image_header;
1472  ECAT63_scanheader scan_header;
1473  ECAT63_attnheader attn_header;
1474  ECAT63_normheader norm_header;
1475  float scale=1.0;
1476  short int *sptr;
1477  char *mdata=NULL, *mptr;
1478  int *iptr;
1479 
1480 
1481  if(IMG_TEST) printf("\nimgReadEcat63Frame(%s, %d, *img, %d)\n",
1482  fname, frame_to_read, frame_index);
1483 
1484  /* Check the input */
1485  if(img==NULL) return STATUS_FAULT;
1486  if(img->status!=IMG_STATUS_OCCUPIED) return STATUS_FAULT;
1487  imgSetStatus(img, STATUS_FAULT);
1488  if(fname==NULL) return STATUS_FAULT;
1489  if(frame_index<0 || frame_index>img->dimt-1) return STATUS_FAULT;
1490  if(frame_to_read<1) return STATUS_FAULT;
1491 
1492  /* Open the file */
1493  if((fp=fopen(fname, "rb")) == NULL) return STATUS_NOFILE;
1494 
1495  /* Read main header */
1496  if((ret=ecat63ReadMainheader(fp, &main_header))) {
1497  fclose(fp); return STATUS_NOMAINHEADER;}
1498 
1499  /* Read matrix list and nr and sort it */
1500  ecat63InitMatlist(&mlist);
1501  ret=ecat63ReadMatlist(fp, &mlist);
1502  if(ret) {fclose(fp); return STATUS_NOMATLIST;}
1503  if(mlist.matrixNr<=0) {
1504  fclose(fp); ecat63EmptyMatlist(&mlist); return STATUS_INVALIDMATLIST;}
1505  /* Make sure that plane and frame numbers are continuous */
1506  ecat63GatherMatlist(&mlist, 1, 1, 1, 1);
1507  ecat63SortMatlistByFrame(&mlist);
1508 
1509  /* Calculate and check the size of one data matrix */
1510  ret=ecat63GetMatrixBlockSize(&mlist, &blkNr);
1511  if(ret) {ecat63EmptyMatlist(&mlist); fclose(fp); return ret;}
1512  /* And allocate memory for the raw data blocks */
1513  if(IMG_TEST>6) printf("allocating memory for %d blocks\n", blkNr);
1514  mdata=(char*)malloc(blkNr*MatBLKSIZE);
1515  if(mdata==NULL) {fclose(fp); ecat63EmptyMatlist(&mlist); return STATUS_NOMEMORY;}
1516 
1517  /* Read all matrices that belong to the required frame */
1518  ret=0; seqplane=-1;
1519  for(m=0; m<mlist.matrixNr; m++) if(mlist.matdir[m].matstat==1) {
1520  /* get frame and plane */
1521  mat_numdoc(mlist.matdir[m].matnum, &matval);
1522  plane=matval.plane;
1523  if(main_header.num_frames>=main_header.num_gates) frame=matval.frame;
1524  else frame=matval.gate; /* printf("frame=%d plane=%d\n", frame, plane); */
1525  if(frame!=frame_to_read) continue; /* do not process other frames */
1526  seqplane++;
1527  if(img->planeNumber[seqplane]<1) {
1528  img->planeNumber[seqplane]=plane;
1529  } else if(img->planeNumber[seqplane]!=plane) {
1530  fclose(fp); ecat63EmptyMatlist(&mlist); free(mdata);
1531  return STATUS_MISSINGMATRIX;
1532  }
1533 
1534  /* Read the subheader */
1535  if(IMG_TEST>4) printf("reading subheader for matrix %d,%d\n", frame, plane);
1536  if(main_header.file_type==IMAGE_DATA) {
1537  ret=ecat63ReadImageheader(fp, mlist.matdir[m].strtblk, &image_header);
1538  scale=image_header.quant_scale;
1539  if(image_header.ecat_calibration_fctr>0.0
1540  && fabs(main_header.calibration_factor/image_header.ecat_calibration_fctr-1.0)>0.0001)
1541  scale*=image_header.ecat_calibration_fctr;
1542  local_data_type=image_header.data_type;
1543  img->start[frame_index]=image_header.frame_start_time/1000.;
1544  img->end[frame_index]=img->start[frame_index]+image_header.frame_duration/1000.;
1545  img->mid[frame_index]=0.5*(img->start[frame_index]+img->end[frame_index]);
1546  if(image_header.decay_corr_fctr>1.0)
1547  img->decayCorrFactor[frame_index]=image_header.decay_corr_fctr;
1548  } else if(main_header.file_type==RAW_DATA) {
1549  ret=ecat63ReadScanheader(fp, mlist.matdir[m].strtblk, &scan_header);
1550  scale=scan_header.scale_factor;
1551  if(scan_header.loss_correction_fctr>0.0) scale*=scan_header.loss_correction_fctr;
1552  local_data_type=scan_header.data_type;
1553  img->start[frame_index]=scan_header.frame_start_time/1000.;
1554  img->end[frame_index]=img->start[frame_index]+scan_header.frame_duration/1000.;
1555  img->mid[frame_index]=0.5*(img->start[frame_index]+img->end[frame_index]);
1556  img->sampleDistance=10.0*scan_header.sample_distance;
1557  img->prompts[frame_index]=(float)scan_header.prompts;
1558  img->randoms[frame_index]=(float)scan_header.delayed;
1559  } else if(main_header.file_type==ATTN_DATA) {
1560  ret=ecat63ReadAttnheader(fp, mlist.matdir[m].strtblk, &attn_header);
1561  scale=attn_header.scale_factor;
1562  local_data_type=attn_header.data_type;
1563  } else if(main_header.file_type==NORM_DATA) {
1564  ret=ecat63ReadNormheader(fp, mlist.matdir[m].strtblk, &norm_header);
1565  scale=norm_header.scale_factor;
1566  local_data_type=norm_header.data_type;
1567  } else
1568  img->_dataType=-1;
1569  if(ret) {
1570  ecat63EmptyMatlist(&mlist); free(mdata); fclose(fp);
1571  return STATUS_NOSUBHEADER;
1572  }
1573  if(main_header.calibration_factor>0.0) scale*=main_header.calibration_factor;
1574  if(IMG_TEST>6) printf("scale=%e datatype=%d\n", scale, img->_dataType);
1575  /* Read the pixel data */
1576  if(IMG_TEST>4) printf("reading matrix data\n");
1577  ret=ecat63ReadMatdata(fp, mlist.matdir[m].strtblk+1,
1578  mlist.matdir[m].endblk-mlist.matdir[m].strtblk,
1579  mdata, local_data_type);
1580  if(ret) {
1581  ecat63EmptyMatlist(&mlist); free(mdata); fclose(fp);
1582  return STATUS_MISSINGMATRIX;
1583  }
1584  if(IMG_TEST>4) printf("converting matrix data to floats\n");
1585  if(local_data_type==BYTE_TYPE) {
1586  for(yi=0, mptr=mdata; yi<img->dimy; yi++)
1587  for(xi=0; xi<img->dimx; xi++, mptr++) {
1588  img->m[seqplane][yi][xi][frame_index]=scale*(float)(*mptr);
1589  }
1590  } else if(local_data_type==VAX_I2 || local_data_type==SUN_I2) {
1591  for(yi=0, mptr=mdata; yi<img->dimy; yi++)
1592  for(xi=0; xi<img->dimx; xi++, mptr+=2) {
1593  sptr=(short int*)mptr;
1594  img->m[seqplane][yi][xi][frame_index]=scale*(float)(*sptr);
1595  }
1596  } else if(local_data_type==VAX_I4 || local_data_type==SUN_I4) {
1597  for(yi=0, mptr=mdata; yi<img->dimy; yi++)
1598  for(xi=0; xi<img->dimx; xi++, mptr+=4) {
1599  iptr=(int*)mptr;
1600  img->m[seqplane][yi][xi][frame_index]=1.0; /* scale*(float)(*iptr); */
1601  }
1602  } else if(local_data_type==VAX_R4 || local_data_type==IEEE_R4) {
1603  for(yi=0, mptr=mdata; yi<img->dimy; yi++)
1604  for(xi=0; xi<img->dimx; xi++, mptr+=4) {
1605  memcpy(&img->m[seqplane][yi][xi][frame_index], mptr, 4);
1606  img->m[seqplane][yi][xi][frame_index]*=scale;
1607  }
1608  }
1609  } /* next matrix */
1610  if(IMG_TEST>3) printf("end of matrices.\n");
1611 
1612  free(mdata);
1613  ecat63EmptyMatlist(&mlist);
1614  fclose(fp);
1615 
1616  /* seqplane is <0 only if this frame did not exist at all; return -1 */
1617  if(IMG_TEST>4) printf("last_seqplane := %d.\n", seqplane);
1618  if(seqplane<0) {imgSetStatus(img, STATUS_NOMATRIX); return STATUS_NOMATRIX;}
1619 
1620  /* check that correct number of planes was read */
1621  if(seqplane+1 != img->dimz) {
1623  return STATUS_MISSINGMATRIX;
1624  }
1625 
1626  /* For saving IMG data, only 2-byte VAX data types are allowed, so change it now */
1627  if(img->_dataType==VAX_I4 || img->_dataType==VAX_R4) img->_dataType=VAX_I2;
1628 
1629  /* All went well */
1630  imgSetStatus(img, STATUS_OK);
1631  return STATUS_OK;
1632 }
1633 /*****************************************************************************/
1634 
1635 /*****************************************************************************/
1656 int imgWriteEcat63Frame(const char *fname, int frame_to_write, IMG *img, int frame_index) {
1657  IMG test_img;
1658  int ret=0, pxlNr, zi, xi, yi, matrixId;
1659  ECAT63_mainheader main_header;
1660  ECAT63_imageheader image_header;
1661  ECAT63_scanheader scan_header;
1662  void *sub_header=NULL;
1663  FILE *fp;
1664  float *fdata=NULL, *fptr;
1665 
1666 
1667  if(IMG_TEST) printf("\nimgWriteEcat63Frame(%s, %d, *img, %d)\n",
1668  fname, frame_to_write, frame_index);
1669  if(IMG_TEST>1) ECAT63_TEST=IMG_TEST-1;
1670 
1671  /*
1672  * Check the input
1673  */
1674  if(fname==NULL) return STATUS_FAULT;
1675  if(img==NULL) return STATUS_FAULT;
1676  if(img->status!=IMG_STATUS_OCCUPIED) return STATUS_FAULT;
1677  if(frame_to_write<0) return STATUS_FAULT;
1678  if(frame_index<0 || frame_index>=img->dimt) return STATUS_FAULT;
1679  if(img->_fileFormat!=IMG_E63) return STATUS_FAULT;
1680 
1681  /* Initiate headers */
1682  memset(&main_header, 0, sizeof(ECAT63_mainheader));
1683  memset(&image_header,0, sizeof(ECAT63_imageheader));
1684  memset(&scan_header, 0, sizeof(ECAT63_scanheader));
1685  imgInit(&test_img);
1686 
1687 
1688  /*
1689  * If file does not exist, then create it with new mainheader,
1690  * and if it does exist, then read and check header information.
1691  * Create or edit mainheader to contain correct frame nr.
1692  * In any case, leave file pointer open for write.
1693  */
1694  if(access(fname, 0) == -1) { /* file does not exist*/
1695 
1696  /* Set main header */
1697  imgSetEcat63MHeader(img, &main_header);
1698  if(IMG_TEST>2) {
1699  ecat63PrintMainheader(&main_header, stdout);
1700  }
1701  if(frame_to_write==0) frame_to_write=1;
1702  main_header.num_frames=frame_to_write;
1703 
1704  /* Open file, write main header and initiate matrix list */
1705  fp=ecat63Create(fname, &main_header); if(fp==NULL) return STATUS_NOWRITEPERM;
1706 
1707  } else { /* file exists*/
1708 
1709  /* Read header information for checking */
1710  ret=imgReadEcat63Header(fname, &test_img); if(ret!=0) return ret;
1711  /* Check that file format is the same */
1712  if(img->_fileFormat!=test_img._fileFormat || img->type!=test_img.type)
1713  return STATUS_WRONGFILETYPE;
1714  /* Check that matrix sizes are the same */
1715  if(img->dimz!=test_img.dimz || img->dimx!=test_img.dimx ||
1716  img->dimy!=test_img.dimy)
1717  return STATUS_VARMATSIZE;
1718  imgEmpty(&test_img);
1719 
1720  /* Open the file for read and write */
1721  if((fp=fopen(fname, "r+b"))==NULL) return STATUS_NOWRITEPERM;
1722 
1723  /* Read the mainheader, set new frame number, and write it back */
1724  if((ret=ecat63ReadMainheader(fp, &main_header))!=0) return STATUS_NOMAINHEADER;
1725  if(frame_to_write==0) frame_to_write=main_header.num_frames+1;
1726  if(main_header.num_frames<frame_to_write)
1727  main_header.num_frames=frame_to_write;
1728  if((ret=ecat63WriteMainheader(fp, &main_header))!=0) return STATUS_NOWRITEPERM;
1729 
1730  }
1731  if(IMG_TEST>2) {
1732  printf("frame_to_write := %d\n", frame_to_write);
1733  }
1734 
1735  /* Allocate memory for matrix float data */
1736  pxlNr=img->dimx*img->dimy*img->dimz;
1737  fdata=(float*)malloc(pxlNr*sizeof(float));
1738  if(fdata==NULL) {fclose(fp); return STATUS_NOMEMORY;}
1739 
1740  /* Set matrix subheader */
1741  if(img->type==IMG_TYPE_RAW) sub_header=(void*)&scan_header;
1742  else if(img->type==IMG_TYPE_IMAGE) sub_header=&image_header;
1743  else {fclose(fp); return STATUS_FAULT;}
1744  imgSetEcat63SHeader(img, sub_header);
1745 
1746  /* Copy matrix pixel values to fdata */
1747  fptr=fdata;
1748  for(zi=0; zi<img->dimz; zi++)
1749  for(yi=0; yi<img->dimy; yi++) for(xi=0; xi<img->dimx; xi++)
1750  *fptr++=img->m[zi][yi][xi][frame_index];
1751 
1752  /* Write subheader and data, and set the rest of subheader contents */
1753  fptr=fdata; ret=-1;
1754  for(zi=0; zi<img->dimz; zi++, fptr+=img->dimx*img->dimy) {
1755  /* Create new matrix id (i.e. matnum) */
1756  matrixId=mat_numcod(frame_to_write, img->planeNumber[zi], 1, 0, 0);
1757  if(img->type==IMG_TYPE_RAW) {
1758  scan_header.frame_start_time=
1759  (int)temp_roundf(1000.*img->start[frame_index]);
1760  scan_header.frame_duration=
1761  (int)temp_roundf(1000.*(img->end[frame_index]-img->start[frame_index]));
1762  scan_header.prompts=temp_roundf(img->prompts[frame_index]);
1763  scan_header.delayed=temp_roundf(img->randoms[frame_index]);
1764  ret=ecat63WriteScanMatrix(fp, matrixId, &scan_header, fptr);
1765  } else {
1766  image_header.frame_start_time=
1767  (int)temp_roundf(1000.*img->start[frame_index]);
1768  image_header.frame_duration=
1769  (int)temp_roundf(1000.*(img->end[frame_index]-img->start[frame_index]));
1770  if(img->decayCorrected)
1771  image_header.decay_corr_fctr=img->decayCorrFactor[frame_index];
1772  ret=ecat63WriteImageMatrix(fp, matrixId, &image_header, fptr);
1773  }
1774  if(ret!=0) break;
1775  } /* next plane*/
1776  free(fdata); fclose(fp);
1777  if(ret) return STATUS_DISKFULL;
1778 
1779  return STATUS_OK;
1780 }
1781 /*****************************************************************************/
1782 
1783 /*****************************************************************************/
1790 void imgSetEcat63SHeader(IMG *img, void *h) {
1791  ECAT63_imageheader *image_header;
1792  ECAT63_scanheader *scan_header;
1793 
1794  if(img->type==IMG_TYPE_RAW) {
1795  scan_header=(ECAT63_scanheader*)h;
1796  scan_header->data_type=VAX_I2;
1797  scan_header->dimension_1=img->dimx;
1798  scan_header->dimension_2=img->dimy;
1799  scan_header->frame_duration_sec=1;
1800  scan_header->scale_factor=1.0;
1801  scan_header->frame_start_time=0;
1802  scan_header->frame_duration=1000;
1803  scan_header->loss_correction_fctr=1.0;
1804  } else {
1805  image_header=(ECAT63_imageheader*)h;
1806  image_header->data_type=VAX_I2;
1807  image_header->num_dimensions=2;
1808  image_header->dimension_1=img->dimx;
1809  image_header->dimension_2=img->dimy;
1810  image_header->recon_scale=img->zoom;
1811  image_header->quant_scale=1.0;
1812  image_header->slice_width=img->sizez/10.;
1813  image_header->pixel_size=img->sizex/10.;
1814  image_header->frame_start_time=0;
1815  image_header->frame_duration=1000;
1816  image_header->plane_eff_corr_fctr=1.0;
1817  image_header->decay_corr_fctr=1.0;
1818  image_header->loss_corr_fctr=1.0;
1819  image_header->ecat_calibration_fctr=1.0;
1820  image_header->well_counter_cal_fctr=1.0;
1821  image_header->quant_units=imgUnitToEcat6(img);
1822  }
1823 }
1824 /*****************************************************************************/
1825 
1826 /*****************************************************************************/
#define SUN_I4
Definition: ecat63.h:36
char ecat63errmsg[128]
Definition: ecat63.h:50
#define ATTN_DATA
Definition: ecat63.h:40
int ECAT63_TEST
Definition: ecat63.h:52
#define VAX_R4
Definition: ecat63.h:33
#define ECAT63_SYSTEM_TYPE_DEFAULT
Definition: ecat63.h:43
#define IEEE_R4
Definition: ecat63.h:34
#define RAW_DATA
Definition: ecat63.h:38
#define BYTE_TYPE
Definition: ecat63.h:30
#define VAX_I4
Definition: ecat63.h:32
#define NORM_DATA
Definition: ecat63.h:41
#define IMAGE_DATA
Definition: ecat63.h:39
#define VAX_I2
Definition: ecat63.h:31
#define MatBLKSIZE
Definition: ecat63.h:27
#define SUN_I2
Definition: ecat63.h:35
int ecat63GetPlaneAndFrameNr(MATRIXLIST *mlist, ECAT63_mainheader *h, int *plane_nr, int *frame_nr)
Definition: ecat63ml.c:414
void ecat63InitMatlist(MATRIXLIST *mlist)
Definition: ecat63ml.c:69
void ecat63EmptyMatlist(MATRIXLIST *mlist)
Definition: ecat63ml.c:80
int ecat63GatherMatlist(MATRIXLIST *ml, short int do_planes, short int do_frames, short int do_gates, short int do_beds)
Definition: ecat63ml.c:519
int ecat63ReadMatlist(FILE *fp, MATRIXLIST *ml)
Definition: ecat63ml.c:97
int mat_numcod(int frame, int plane, int gate, int data, int bed)
Definition: ecat63ml.c:266
int ecat63GetMatrixBlockSize(MATRIXLIST *mlist, int *blk_nr)
Definition: ecat63ml.c:382
int ecat63DeleteLateFrames(MATRIXLIST *ml, int frame_nr)
Definition: ecat63ml.c:360
void ecat63SortMatlistByPlane(MATRIXLIST *ml)
Definition: ecat63ml.c:291
void ecat63PrintMatlist(MATRIXLIST *ml)
Definition: ecat63ml.c:160
void mat_numdoc(int matnum, Matval *matval)
Definition: ecat63ml.c:276
void ecat63SortMatlistByFrame(MATRIXLIST *ml)
Definition: ecat63ml.c:316
void ecat63PrintMainheader(ECAT63_mainheader *h, FILE *fp)
Definition: ecat63p.c:62
void ecat63PrintImageheader(ECAT63_imageheader *h, FILE *fp)
Definition: ecat63p.c:115
void ecat63PrintScanheader(ECAT63_scanheader *h, FILE *fp)
Definition: ecat63p.c:152
int ecat63ReadNormheader(FILE *fp, int blk, ECAT63_normheader *h)
Definition: ecat63r.c:375
int ecat63ReadMatdata(FILE *fp, int strtblk, int blkNr, char *data, int dtype)
Definition: ecat63r.c:432
int ecat63ReadScanheader(FILE *fp, int blk, ECAT63_scanheader *h)
Definition: ecat63r.c:296
int ecat63ReadAttnheader(FILE *fp, int blk, ECAT63_attnheader *h)
Definition: ecat63r.c:238
int ecat63ReadImageheader(FILE *fp, int blk, ECAT63_imageheader *h)
Definition: ecat63r.c:152
int ecat63ReadMainheader(FILE *fp, ECAT63_mainheader *h)
Definition: ecat63r.c:50
int ecat63WriteScan(FILE *fp, int matnum, ECAT63_scanheader *h, void *data)
Definition: ecat63w.c:478
int ecat63WriteScanMatrix(FILE *fp, int matnum, ECAT63_scanheader *h, float *fdata)
Definition: ecat63w.c:784
int ecat63WriteImageMatrix(FILE *fp, int matnum, ECAT63_imageheader *h, float *fdata)
Definition: ecat63w.c:700
FILE * ecat63Create(const char *fname, ECAT63_mainheader *h)
Definition: ecat63w.c:386
int ecat63WriteImage(FILE *fp, int matnum, ECAT63_imageheader *h, void *data)
Definition: ecat63w.c:429
int ecat63WriteMainheader(FILE *fp, ECAT63_mainheader *h)
Definition: ecat63w.c:73
void imgInfo(IMG *image)
Definition: img.c:414
int imgAllocate(IMG *image, int planes, int rows, int columns, int frames)
Definition: img.c:285
void imgSetStatus(IMG *img, int status_index)
Definition: img.c:399
void imgEmpty(IMG *image)
Definition: img.c:216
void imgInit(IMG *image)
Definition: img.c:163
#define IMG_TYPE_RAW
Definition: img.h:81
@ STATUS_NOMATLIST
Definition: img.h:120
@ STATUS_OK
Definition: img.h:118
@ STATUS_DISKFULL
Definition: img.h:119
@ STATUS_WRONGFILETYPE
Definition: img.h:124
@ STATUS_NOWRITEPERM
Definition: img.h:119
@ STATUS_NOFILE
Definition: img.h:118
@ STATUS_NOMATRIX
Definition: img.h:121
@ STATUS_MISSINGMATRIX
Definition: img.h:119
@ STATUS_INVALIDMATLIST
Definition: img.h:120
@ STATUS_NOSUBHEADER
Definition: img.h:121
@ STATUS_VARMATSIZE
Definition: img.h:120
@ STATUS_UNSUPPORTED
Definition: img.h:119
@ STATUS_FAULT
Definition: img.h:118
@ STATUS_NOMAINHEADER
Definition: img.h:120
@ STATUS_NOMEMORY
Definition: img.h:118
int IMG_TEST
Definition: img.h:128
#define IMG_STATUS_OCCUPIED
Definition: img.h:73
#define IMG_UNKNOWN
Definition: img.h:84
#define IMG_STATUS_INITIALIZED
Definition: img.h:72
#define IMG_E63
Definition: img.h:85
#define IMG_STATUS_UNINITIALIZED
Definition: img.h:71
#define IMG_TYPE_IMAGE
Definition: img.h:80
int ecat63ReadPlaneToImg(const char *fname, IMG *img)
Definition: img_e63.c:568
int imgGetEcat63Fileformat(ECAT63_mainheader *h)
Definition: img_e63.c:1260
int imgEcat63Supported(ECAT63_mainheader *h)
Definition: img_e63.c:1140
int ecat63WriteAllImg(const char *fname, IMG *img)
Definition: img_e63.c:374
void imgSetEcat63SHeader(IMG *img, void *h)
Definition: img_e63.c:1790
int imgWriteEcat63Frame(const char *fname, int frame_to_write, IMG *img, int frame_index)
Definition: img_e63.c:1656
int imgReadEcat63FirstFrame(const char *fname, IMG *img)
Definition: img_e63.c:1422
int ecat63AddImg(const char *fname, IMG *img)
Definition: img_e63.c:886
void imgGetEcat63MHeader(IMG *img, ECAT63_mainheader *h)
Definition: img_e63.c:1157
int imgReadEcat63Header(const char *fname, IMG *img)
Definition: img_e63.c:1289
void imgSetEcat63MHeader(IMG *img, ECAT63_mainheader *h)
Definition: img_e63.c:1208
int imgReadEcat63Frame(const char *fname, int frame_to_read, IMG *img, int frame_index)
Definition: img_e63.c:1464
int ecat63ReadAllToImg(const char *fname, IMG *img)
Definition: img_e63.c:77
char * imgIsotope(IMG *img)
Definition: imgdecay.c:110
int imgUnitToEcat6(IMG *img)
Definition: imgunit.c:233
void imgUnitFromEcat(IMG *img, int ecat_unit)
Definition: imgunit.c:160
Definition: img.h:156
float sizex
Definition: img.h:208
unsigned short int dimx
Definition: img.h:261
char type
Definition: img.h:198
char patientName[32]
Definition: img.h:176
char studyDescription[32]
Definition: img.h:192
float sampleDistance
Definition: img.h:206
float **** m
Definition: img.h:274
float transaxialFOV
Definition: img.h:204
char unit
Definition: img.h:172
char status
Definition: img.h:164
time_t scanStart
Definition: img.h:186
int _fileFormat
Definition: img.h:229
char decayCorrected
Definition: img.h:184
float * prompts
Definition: img.h:306
char userProcessCode[11]
Definition: img.h:190
unsigned short int dimt
Definition: img.h:259
int _dataType
Definition: img.h:226
int * planeNumber
Definition: img.h:284
char patientID[16]
Definition: img.h:178
int scanner
Definition: img.h:231
float sizey
Definition: img.h:210
float * start
Definition: img.h:290
unsigned short int dimz
Definition: img.h:265
unsigned short int dimy
Definition: img.h:263
float * end
Definition: img.h:292
float zoom
Definition: img.h:200
char radiopharmaceutical[32]
Definition: img.h:180
float * decayCorrFactor
Definition: img.h:314
float isotopeHalflife
Definition: img.h:182
char studyNr[MAX_STUDYNR_LEN+1]
Definition: img.h:174
float * randoms
Definition: img.h:308
float axialFOV
Definition: img.h:202
float * mid
Definition: img.h:294
float sizez
Definition: img.h:212
MatDir * matdir
Definition: ecat63.h:65
int matrixNr
Definition: ecat63.h:63
Definition: ecat63.h:55
int matstat
Definition: ecat63.h:59
int endblk
Definition: ecat63.h:58
int matnum
Definition: ecat63.h:56
int strtblk
Definition: ecat63.h:57
Definition: ecat63.h:68
int frame
Definition: ecat63.h:69
int gate
Definition: ecat63.h:69
int plane
Definition: ecat63.h:69
float scale_factor
Definition: ecat63.h:157
short int data_type
Definition: ecat63.h:155
short int dimension_1
Definition: ecat63.h:156
short int dimension_2
Definition: ecat63.h:156
short int data_type
Definition: ecat63.h:107
int frame_start_time
Definition: ecat63.h:111
short int quant_units
Definition: ecat63.h:119
short int num_dimensions
Definition: ecat63.h:107
float decay_corr_fctr
Definition: ecat63.h:118
float well_counter_cal_fctr
Definition: ecat63.h:121
short int image_max
Definition: ecat63.h:109
float pixel_size
Definition: ecat63.h:110
float plane_eff_corr_fctr
Definition: ecat63.h:117
short int dimension_1
Definition: ecat63.h:107
float recon_scale
Definition: ecat63.h:108
float ecat_calibration_fctr
Definition: ecat63.h:121
float loss_corr_fctr
Definition: ecat63.h:118
float slice_width
Definition: ecat63.h:110
float quant_scale
Definition: ecat63.h:108
short int dimension_2
Definition: ecat63.h:107
short int image_min
Definition: ecat63.h:109
short int num_frames
Definition: ecat63.h:97
short int scan_start_minute
Definition: ecat63.h:81
short int scan_start_hour
Definition: ecat63.h:81
short int scan_start_day
Definition: ecat63.h:80
float transaxial_fov
Definition: ecat63.h:87
short int sw_version
Definition: ecat63.h:75
short int num_gates
Definition: ecat63.h:97
char radiopharmaceutical[32]
Definition: ecat63.h:84
char study_description[32]
Definition: ecat63.h:94
short int data_type
Definition: ecat63.h:76
char patient_name[32]
Definition: ecat63.h:91
short int system_type
Definition: ecat63.h:77
float isotope_halflife
Definition: ecat63.h:83
float calibration_factor
Definition: ecat63.h:89
char patient_id[16]
Definition: ecat63.h:91
short int scan_start_year
Definition: ecat63.h:80
short int scan_start_month
Definition: ecat63.h:80
float axial_fov
Definition: ecat63.h:87
char isotope_code[8]
Definition: ecat63.h:82
short int calibration_units
Definition: ecat63.h:90
short int num_bed_pos
Definition: ecat63.h:97
char study_name[12]
Definition: ecat63.h:91
short int scan_start_second
Definition: ecat63.h:81
char user_process_code[10]
Definition: ecat63.h:101
short int file_type
Definition: ecat63.h:78
short int num_planes
Definition: ecat63.h:97
float plane_separation
Definition: ecat63.h:98
short int dimension_1
Definition: ecat63.h:148
float scale_factor
Definition: ecat63.h:149
short int dimension_2
Definition: ecat63.h:148
short int data_type
Definition: ecat63.h:147
int frame_duration
Definition: ecat63.h:141
float loss_correction_fctr
Definition: ecat63.h:142
short int dimension_2
Definition: ecat63.h:129
short int scan_max
Definition: ecat63.h:136
int frame_start_time
Definition: ecat63.h:141
short int frame_duration_sec
Definition: ecat63.h:133
short int scan_min
Definition: ecat63.h:136
float sample_distance
Definition: ecat63.h:131
float scale_factor
Definition: ecat63.h:135
short int data_type
Definition: ecat63.h:128
short int dimension_1
Definition: ecat63.h:129