img_ana.c
Go to the documentation of this file.
1 /******************************************************************************
2 
3  Copyright (c) 2007,2008,2009 Turku PET Centre
4 
5  Library: img_ana.c
6  Description: I/O routines for IMG data from/to Analyze 7.5 format.
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  Prompts and randoms are read from SIF data to IMG when reading Analyze.
29  2007-02-27 VO
30  Bug correction in imgWriteAnalyze(): fp was not closed in all errors.
31  2007-03-25 VO
32  Added functions imgReadAnalyzeHeader(), imgGetAnalyzeHeader(),
33  imgSetAnalyzeHeader(), imgReadAnalyzeFrame(), imgReadAnalyzeFirstFrame(),
34  and imgWriteAnalyzeFrame().
35  2007-17-07 Harri Merisaari
36  Modified for optional ANSi compatibility
37  2007-09-10 VO
38  Return value of localtime() is checked.
39  2008-07-07 VO
40  If information on decay correction is not found in Analyze header, then
41  it is assumed that image data is corrected for decay; previously assumed
42  that image was NOT corrected for decay.
43  2009-12-10 VO
44  strcpy() replaced with strncpy() in filling header fields.
45 
46 
47 ******************************************************************************/
48 #include <stdio.h>
49 #include <stdlib.h>
50 #include <unistd.h>
51 #include <math.h>
52 #include <string.h>
53 #include <time.h>
54 /*****************************************************************************/
55 #include "petc99.h"
56 #include "swap.h"
57 #include "halflife.h"
58 /*****************************************************************************/
59 #include "include/img.h"
60 #include "include/analyze.h"
61 #include "include/imgmax.h"
62 #include "include/imgdecay.h"
63 #include "include/sif.h"
64 #include "include/imgfile.h"
65 /*****************************************************************************/
66 
67 /*****************************************************************************/
83 int imgReadAnalyze(const char *dbname, IMG *img) {
84  FILE *fp;
85  int ret, fi, pi, xi, yi;
86  float *fdata=NULL, *fptr;
87  ANALYZE_DSR dsr;
88  char datfile[FILENAME_MAX], hdrfile[FILENAME_MAX], siffile[FILENAME_MAX];
89  char buf[128], *cptr;
90  int dimNr, dimx, dimy, dimz=1, dimt=1, pxlNr=0;
91  SIF sif;
92  struct tm *st;
93 
94 
95  if(IMG_TEST) printf("imgReadAnalyze(%s, *img)\n", dbname);
96 
97  /* Check the arguments */
98  imgSetStatus(img, STATUS_OK);
99  if(img==NULL || img->status!=IMG_STATUS_INITIALIZED) {
100  imgSetStatus(img, STATUS_FAULT); return(2);}
101  if(dbname==NULL || !dbname[0]) {imgSetStatus(img, STATUS_FAULT); return(1);}
102 
103  /* Make the image and header filenames */
104  if(anaExists(dbname)==0) {
105  /* Check if filename was given accidentally with extension */
106  strcpy(datfile, dbname); cptr=strrchr(datfile, '.');
107  if(cptr!=NULL && (strcmp(cptr, ".img")==0 || strcmp(cptr, ".hdr")==0)) {
108  *cptr=(char)0; strcpy(hdrfile, datfile);
109  if(anaExists(datfile)==0) { /* still not found */
110  imgSetStatus(img, STATUS_NOHEADERFILE); return(3);}
111  strcat(datfile, ".img"); strcat(hdrfile, ".hdr");
112  } else {
113  imgSetStatus(img, STATUS_NOHEADERFILE); return(3);
114  }
115  } else {
116  /* Database name was given and img and hdr files were found */
117  strcpy(datfile, dbname); strcat(datfile, ".img");
118  strcpy(hdrfile, dbname); strcat(hdrfile, ".hdr");
119  }
120 
121  /* Read Analyze header file */
122  ret=anaReadHeader(hdrfile, &dsr);
123  if(ret) {
124  if(ret==1) imgSetStatus(img, STATUS_FAULT);
125  else if(ret==2) imgSetStatus(img, STATUS_NOHEADERFILE);
126  else imgSetStatus(img, STATUS_UNSUPPORTED);
127  return(3);
128  }
129  if(IMG_TEST) anaPrintHeader(&dsr, stdout);
130 
131  /* Open image datafile */
132  if(IMG_TEST) fprintf(stdout, "reading image data %s\n", datfile);
133  if((fp=fopen(datfile, "rb")) == NULL) {imgSetStatus(img, STATUS_NOIMGDATA); return(5);}
134 
135  /* Prepare IMG for Analyze image */
136  /* Get the image dimensions from header */
137  dimNr=dsr.dime.dim[0];
138  if(dimNr<2) {fclose(fp); imgSetStatus(img, STATUS_INVALIDHEADER); return(4);}
139  dimx=dsr.dime.dim[1]; dimy=dsr.dime.dim[2];
140  if(dimNr>2) {dimz=dsr.dime.dim[3]; if(dimNr>3) dimt=dsr.dime.dim[4];}
141  pxlNr=dimx*dimy*dimz;
142  if(pxlNr<1) {fclose(fp); imgSetStatus(img, STATUS_INVALIDHEADER); return(4);}
143  /* Allocate memory for IMG */
144  ret=imgAllocate(img, dimz, dimy, dimx, dimt);
145  if(ret) {fclose(fp); imgSetStatus(img, STATUS_NOMEMORY); return(11);}
146  /* Copy information from Analyze header */
147  img->type=IMG_TYPE_IMAGE;
148  strncpy(img->studyNr, dsr.hist.patient_id, MAX_STUDYNR_LEN);
149  strcpy(img->patientName, dsr.hist.patient_id);
150  img->sizex=dsr.dime.pixdim[1];
151  img->sizey=dsr.dime.pixdim[2];
152  img->sizez=dsr.dime.pixdim[3];
153  /*if(dsr.dime.funused2>1.E-5) img->zoom=dsr.dime.funused2;*/
154  if(dsr.dime.funused3>1.E-5) img->isotopeHalflife=dsr.dime.funused3;
155  for(pi=0; pi<dimz; pi++) img->planeNumber[pi]=pi+1;
156  if(dsr.little) img->_fileFormat=IMG_ANA_L; else img->_fileFormat=IMG_ANA;
157  /* Decay correction */
158  if(strstr(dsr.hist.descrip, "Decay corrected.")!=NULL)
159  img->decayCorrected=1;
160  else if(strstr(dsr.hist.descrip, "No decay correction.")!=NULL)
161  img->decayCorrected=0;
162  else
163  img->decayCorrected=1;
164 
165  /* Allocate memory for one image frame */
166  fdata=malloc(pxlNr*sizeof(float));
167  if(fdata==NULL) {fclose(fp); imgSetStatus(img, STATUS_NOMEMORY); return(12);}
168 
169  /* Read one image frame at a time */
170  for(fi=0; fi<dimt; fi++) {
171  fptr=fdata;
172  ret=anaReadImagedata(fp, &dsr, fi+1, fptr);
173  if(ret) {free(fdata); fclose(fp); imgSetStatus(img, STATUS_NOIMGDATA); return(7);}
174  /* Copy pixel values to IMG */
175  fptr=fdata;
176  if(anaFlipping()==0) { /* no flipping in z-direction */
177  for(pi=0; pi<img->dimz; pi++)
178  for(yi=dimy-1; yi>=0; yi--)
179  for(xi=dimx-1; xi>=0; xi--)
180  img->m[pi][yi][xi][fi]=*fptr++;
181  } else {
182  for(pi=dimz-1; pi>=0; pi--)
183  for(yi=dimy-1; yi>=0; yi--)
184  for(xi=dimx-1; xi>=0; xi--)
185  img->m[pi][yi][xi][fi]=*fptr++;
186  }
187  } /* next frame */
188  free(fdata);
189  fclose(fp);
190 
191  /* Try to read frame time information from SIF file */
192  /* Make filename from database or image data file */
193  strcpy(siffile, dbname); strcat(siffile, ".sif");
194  if(access(siffile, 0) == -1) {
195  strcpy(siffile, datfile); strcat(siffile, ".sif");
196  }
197  /* Check if SIF file is found */
198  if(IMG_TEST) printf("reading SIF file %s\n", siffile);
199  if(access(siffile, 0) == -1) {
200  if(IMG_TEST) printf(" No SIF file; therefore unknown frame times.\n");
201  return(0);
202  }
203  /* If found, then read it */
204  sifInit(&sif); ret=sifRead(siffile, &sif);
205  if(ret) {imgSetStatus(img, STATUS_NOSIFDATA); return(21);}
206  /* Copy scan start time */
207  if(sif.scantime>0) {
208  st=localtime(&sif.scantime);
209  if(st!=NULL) strftime(buf, 128, "%Y-%m-%d %H:%M:%S", st);
210  else strcpy(buf, "1900-01-01 00:00:00");
211  img->scanStart=sif.scantime;
212  }
213  /* Copy frame times */
214  if(sif.frameNr!=img->dimt) {
215  imgSetStatus(img, STATUS_WRONGSIFDATA); sifEmpty(&sif); return(22);}
216  for(fi=0; fi<sif.frameNr; fi++) {
217  img->start[fi]=sif.x1[fi]; img->end[fi]=sif.x2[fi];
218  img->mid[fi]=0.5*(img->start[fi]+img->end[fi]);
219  }
220  /* Copy prompts and randoms */
221  for(fi=0; fi<sif.frameNr; fi++) {
222  img->prompts[fi]=sif.prompts[fi];
223  img->randoms[fi]=sif.randoms[fi];
224  }
225 
226  /* Set isotopeHalflife, if isotope is found */
227  if(strlen(sif.isotope_name)>1) {
228  img->isotopeHalflife=60.0*hlFromIsotope(sif.isotope_name);
229  }
230  sifEmpty(&sif);
231 
232  return(0);
233 }
234 /*****************************************************************************/
235 
236 /*****************************************************************************/
253 int imgWriteAnalyze(const char *dbname, IMG *img) {
254  FILE *fp;
255  int ret, fi, pi, xi, yi, little;
256  float g;
257  ANALYZE_DSR dsr;
258  char datfile[FILENAME_MAX], hdrfile[FILENAME_MAX];
259  const char *cptr;
260  int pxlNr=0;
261  struct tm *st;
262  short int *sdata, *sptr, smin, smax;
263 
264 
265  if(IMG_TEST) printf("imgWriteAnalyze(%s, *img)\n", dbname);
266 
267  /* Check the arguments */
268  imgSetStatus(img, STATUS_OK);
269  if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) {
270  imgSetStatus(img, STATUS_FAULT); return(2);}
271  if(dbname==NULL || !dbname[0]) {imgSetStatus(img, STATUS_FAULT); return(1);}
272 
273  /* Make the image and header filenames */
274  strcpy(datfile, dbname); strcat(datfile, ".img");
275  strcpy(hdrfile, dbname); strcat(hdrfile, ".hdr");
276 
277 
278  /*
279  * Fill Analyze header
280  */
281  if(img->_fileFormat==IMG_ANA_L) dsr.little=1; else dsr.little=0;
282  /* Header key */
283  memset(&dsr.hk, 0, sizeof(ANALYZE_HEADER_KEY));
284  memset(&dsr.dime, 0, sizeof(ANALYZE_HEADER_IMGDIM));
285  memset(&dsr.hist, 0, sizeof(ANALYZE_HEADER_HISTORY));
286  dsr.hk.sizeof_hdr=348;
287  strcpy(dsr.hk.data_type, "");
288  cptr=strrchr(dbname, '/'); if(cptr==NULL) cptr=strrchr(dbname, '\\');
289  if(cptr!=NULL) cptr++; if(cptr==NULL) cptr=dbname;
290  strncpy(dsr.hk.db_name, cptr, 17);
291  dsr.hk.extents=16384;
292  dsr.hk.regular='r';
293  /* Image dimension */
294  dsr.dime.dim[0]=4;
295  dsr.dime.dim[1]=img->dimx;
296  dsr.dime.dim[2]=img->dimy;
297  dsr.dime.dim[3]=img->dimz;
298  dsr.dime.dim[4]=img->dimt;
300  dsr.dime.bitpix=16;
301  dsr.dime.pixdim[0]=0.0;
302  dsr.dime.pixdim[1]=img->sizex;
303  dsr.dime.pixdim[2]=img->sizey;
304  dsr.dime.pixdim[3]=img->sizez;
305  dsr.dime.pixdim[4]=0.0;
306  dsr.dime.funused1=0.0; /* Scale factor is set later */
307  /* dsr.dime.funused2=img->zoom; */ /* Reconstruction zoom */
308  dsr.dime.funused3=img->isotopeHalflife;
309  /* Data history */
310  if(img->decayCorrected==1) strcpy(dsr.hist.descrip, "Decay corrected.");
311  else strcpy(dsr.hist.descrip, "No decay correction.");
312  strncpy(dsr.hist.scannum, img->studyNr, 10);
313  st=localtime(&img->scanStart);
314  if(st!=NULL) {
315  strftime(dsr.hist.exp_date, 10, "%Y-%m-%d", st);
316  strftime(dsr.hist.exp_time, 10, "%H:%M:%S", st);
317  } else {
318  strncpy(dsr.hist.exp_date, "1900-01-01", 10);
319  strncpy(dsr.hist.exp_time, "00:00:00", 10);
320  }
321 
322  /*
323  * Scale data to short int range
324  * Determine and set scale factor and cal_min & cal_max
325  */
326  if(IMG_TEST) printf("scaling data to short ints\n");
327  ret=imgMinMax(img, &dsr.dime.cal_min, &dsr.dime.cal_max);
328  if(ret) {imgSetStatus(img, STATUS_FAULT); return(3);}
329  if(fabs(dsr.dime.cal_min)>fabs(dsr.dime.cal_max)) g=fabs(dsr.dime.cal_min);
330  else g=fabs(dsr.dime.cal_max);
331  if(g<1E-20) g=1.0; else g=32767./g; dsr.dime.funused1=1.0/g;
332  if(IMG_TEST) printf("min=%g max=%g scale_factor=%g\n",
333  dsr.dime.cal_min, dsr.dime.cal_max, dsr.dime.funused1);
334 
335  /* Allocate memory for short int array */
336  pxlNr=(img->dimx)*(img->dimy)*(img->dimz);
337  sdata=malloc(pxlNr*sizeof(short int));
338  if(sdata==NULL) {
340  return 12;
341  }
342 
343  /* Open image data file for write */
344  if((fp=fopen(datfile, "wb")) == NULL) {
346  free(sdata);
347  return 14;
348  }
349 
350  /* Copy and write image matrix data to short int array */
351  /* Data is written one frame at a time */
352  smin=smax=temp_roundf(g*img->m[0][0][0][0]);
353  for(fi=0; fi<img->dimt; fi++) {
354  sptr=sdata;
355  if(anaFlipping()==0) {
356  for(pi=0; pi<img->dimz; pi++)
357  for(yi=img->dimy-1; yi>=0; yi--)
358  for(xi=img->dimx-1; xi>=0; xi--) {
359  *sptr=temp_roundf(g*img->m[pi][yi][xi][fi]);
360  if(*sptr>smax) smax=*sptr; else if(*sptr<smin) smin=*sptr;
361  sptr++;
362  }
363  } else {
364  for(pi=img->dimz-1; pi>=0; pi--)
365  for(yi=img->dimy-1; yi>=0; yi--)
366  for(xi=img->dimx-1; xi>=0; xi--) {
367  *sptr=temp_roundf(g*img->m[pi][yi][xi][fi]);
368  if(*sptr>smax) smax=*sptr; else if(*sptr<smin) smin=*sptr;
369  sptr++;
370  }
371  }
372  /* Change byte order if necessary */
373  little=little_endian();
374  if(little!=dsr.little)
375  swabip(sdata, pxlNr*sizeof(short int));
376  /* Write image data */
377  if(fwrite(sdata, 2, pxlNr, fp) != pxlNr) {
379  free(sdata); fclose(fp);
380  return 15;
381  }
382  }
383  /* Done writing */
384  fclose(fp);
385  free(sdata);
386 
387  if(IMG_TEST) printf("smin=%d smax=%d\n", smin, smax);
388 
389  /* Set header glmin & glmax */
390  dsr.dime.glmin=smin; dsr.dime.glmax=smax;
391 
392  /* Write Analyze header */
393  ret=anaWriteHeader(hdrfile, &dsr);
394  if(ret) {
396  return 21;
397  }
398  imgSetStatus(img, STATUS_OK);
399  return 0;
400 }
401 /*****************************************************************************/
402 
403 /*****************************************************************************/
414 int imgReadAnalyzeHeader(const char *dbname, IMG *img) {
415  char hdrfile[FILENAME_MAX], siffile[FILENAME_MAX];
416  ANALYZE_DSR ana_header;
417  SIF sif;
418  double f;
419  int ret;
420 
421  if(IMG_TEST) printf("\nimgReadAnalyzeHeader(%s, *img)\n", dbname);
422 
423  /* Check the input */
424  if(img==NULL) return STATUS_FAULT;
425  if(img->status!=IMG_STATUS_INITIALIZED) return STATUS_FAULT;
427  if(dbname==NULL) return STATUS_FAULT;
428 
429  /* Determine the names of hdr and sif files */
430  ret=anaDatabaseExists(dbname, hdrfile, NULL, siffile);
431  if(ret==0) return STATUS_NOFILE;
432 
433  /* Read Analyze header file */
434  ret=anaReadHeader(hdrfile, &ana_header);
435  if(ret!=0) {
436  if(IMG_TEST>1) printf("anaReadHeader() return value := %d\n", ret);
437  if(ret==1) return STATUS_FAULT;
438  else if(ret==2) return STATUS_NOHEADERFILE;
439  else return STATUS_UNSUPPORTED;
440  return(STATUS_FAULT);
441  }
442  /* and set IMG contents */
443  ret=imgGetAnalyzeHeader(img, &ana_header);
444  if(ret!=0) {
445  imgSetStatus(img, ret);
446  return(ret);
447  }
448 
449  /* If SIF does not exist, then that's it */
450  if(!siffile[0]) {
451  imgSetStatus(img, STATUS_OK);
452  return STATUS_OK;
453  }
454 
455  /* SIF is available, so read that too */
456  sifInit(&sif); ret=0;
457  if(sifRead(siffile, &sif)!=0) return STATUS_OK;
458  /* Copy scan time */
459  img->scanStart=sif.scantime;
460  /* Study number, if not yet defined */
461  if(!img->studyNr[0] && strlen(sif.studynr)>1 )
462  strncpy(img->studyNr, sif.studynr, MAX_STUDYNR_LEN);
463  /* Isotope half-life, if not yet defined */
464  f=hlFromIsotope(sif.isotope_name);
465  if(img->isotopeHalflife<=0.0 && f>0.0) img->isotopeHalflife=60.0*f;
466  sifEmpty(&sif);
467 
468  return STATUS_OK;
469 }
470 /*****************************************************************************/
471 
472 /*****************************************************************************/
482  int dimNr, dimx, dimy, dimz=1, dimt=1, pxlNr=0;
483 
484  if(IMG_TEST) printf("\nimgGetAnalyzeHeader(*img, *dsr)\n");
485 
486  /* Check the input */
487 
488  if(img==NULL) return STATUS_FAULT;
490  return STATUS_FAULT;
492  if(h==NULL) return STATUS_FAULT;
493 
495 
496  /* Get the image dimensions from header */
497  dimNr=h->dime.dim[0];
498  if(dimNr<2) return STATUS_INVALIDHEADER;
499  dimx=h->dime.dim[1]; dimy=h->dime.dim[2];
500  if(dimNr>2) {dimz=h->dime.dim[3]; if(dimNr>3) dimt=h->dime.dim[4];}
501  pxlNr=dimx*dimy*dimz;
502  if(pxlNr<1) return STATUS_INVALIDHEADER;
503  img->dimx=dimx; img->dimy=dimy; img->dimz=dimz; img->dimt=dimt;
504 
505  /* Copy information from Analyze header */
506  img->type=IMG_TYPE_IMAGE;
507  strncpy(img->studyNr, h->hist.patient_id, MAX_STUDYNR_LEN);
508  strcpy(img->patientName, h->hist.patient_id);
509  img->sizex=h->dime.pixdim[1];
510  img->sizey=h->dime.pixdim[2];
511  img->sizez=h->dime.pixdim[3];
512  /*if(h->dime.funused2>1.E-5) img->zoom=h->dime.funused2;*/
513  if(h->dime.funused3>1.E-5) img->isotopeHalflife=h->dime.funused3;
514  if(h->little) img->_fileFormat=IMG_ANA_L; else img->_fileFormat=IMG_ANA;
515 
516  /* Decay correction */
517  if(strstr(h->hist.descrip, "Decay corrected.")!=NULL)
518  img->decayCorrected=1;
519  else if(strstr(h->hist.descrip, "No decay correction.")!=NULL)
520  img->decayCorrected=0;
521  else
522  img->decayCorrected=1;
523 
524  imgSetStatus(img, STATUS_OK);
525  return STATUS_OK;
526 }
527 /*****************************************************************************/
528 
529 /*****************************************************************************/
542 int imgSetAnalyzeHeader(IMG *img, const char *dbname, ANALYZE_DSR *dsr, float fmin, float fmax) {
543  struct tm *st;
544  char *cptr;
545  float g;
546 
547  if(IMG_TEST) printf("\nimgSetAnalyzeHeader(*img, *dsr)\n");
548 
549  /* Check the input */
550  if(img==NULL) return STATUS_FAULT;
552  return STATUS_FAULT;
554  if(dsr==NULL) return STATUS_FAULT;
555 
556  /* Byte order */
557  if(img->_fileFormat==IMG_ANA_L) dsr->little=1; else dsr->little=0;
558  /* Header key */
559  memset(&dsr->hk, 0, sizeof(ANALYZE_HEADER_KEY));
560  memset(&dsr->dime, 0, sizeof(ANALYZE_HEADER_IMGDIM));
561  memset(&dsr->hist, 0, sizeof(ANALYZE_HEADER_HISTORY));
562  dsr->hk.sizeof_hdr=348;
563  strcpy(dsr->hk.data_type, "");
564  cptr=strrchr(dbname, '/'); if(cptr==NULL) cptr=strrchr(dbname, '\\');
565  if(cptr!=NULL) cptr++; if(cptr==NULL) cptr=(char*)dbname;
566  strncpy(dsr->hk.db_name, cptr, 17);
567  dsr->hk.extents=16384;
568  dsr->hk.regular='r';
569  /* Image dimension */
570  dsr->dime.dim[0]=4;
571  dsr->dime.dim[1]=img->dimx;
572  dsr->dime.dim[2]=img->dimy;
573  dsr->dime.dim[3]=img->dimz;
574  dsr->dime.dim[4]=img->dimt;
576  dsr->dime.bitpix=16;
577  dsr->dime.pixdim[0]=0.0;
578  dsr->dime.pixdim[1]=img->sizex;
579  dsr->dime.pixdim[2]=img->sizey;
580  dsr->dime.pixdim[3]=img->sizez;
581  dsr->dime.pixdim[4]=0.0;
582  dsr->dime.funused1=0.0; /* Scale factor is set later */
583  /* dsr.dime.funused2=img->zoom; */ /* Reconstruction zoom */
584  dsr->dime.funused3=img->isotopeHalflife;
585  /* Data history */
586  if(img->decayCorrected==1) strcpy(dsr->hist.descrip, "Decay corrected.");
587  else strcpy(dsr->hist.descrip, "No decay correction.");
588  strncpy(dsr->hist.scannum, img->studyNr, 10);
589  st=localtime(&img->scanStart);
590  if(st!=NULL) {
591  strftime(dsr->hist.exp_date, 10, "%Y-%m-%d", st);
592  strftime(dsr->hist.exp_time, 10, "%H:%M:%S", st);
593  } else {
594  strcpy(dsr->hist.exp_date, "1900-01-01");
595  strcpy(dsr->hist.exp_time, "00:00:00");
596  }
597  /* Determine and set scale factor and cal_min & cal_max */
598  if(fmin<fmax) {
599  dsr->dime.cal_min=fmin; dsr->dime.cal_max=fmax;
600  } else { /* not given in function call, try to find those here */
601  if(img->status==IMG_STATUS_OCCUPIED &&
602  imgMinMax(img, &dsr->dime.cal_min, &dsr->dime.cal_max)==0) {}
603  else return STATUS_FAULT;
604  }
605  if(fabs(dsr->dime.cal_min) > fabs(dsr->dime.cal_max)) g = fabs(dsr->dime.cal_min);
606  else g = fabs(dsr->dime.cal_max);
607  /* if(fabs(dsr->dime.cal_min)>fabs(dsr->dime.cal_max)) g=fabs(dsr->dime.cal_min); */
608  /* else g=fabs(dsr->dime.cal_max); */
609  if(g<1E-20) g=1.0; else g=32767./g; dsr->dime.funused1=1.0/g;
610  /* Set header glmin & glmax */
611  dsr->dime.glmin=temp_roundf(fmin*g); dsr->dime.glmax=temp_roundf(fmax*g);
612  /* printf("glmin=%d\n", dsr->dime.glmin); */
613  /* printf("glmax=%d\n", dsr->dime.glmax); */
614 
615  imgSetStatus(img, STATUS_OK);
616  return STATUS_OK;
617 }
618 /*****************************************************************************/
619 
620 /*****************************************************************************/
629 int imgReadAnalyzeFirstFrame(const char *fname, IMG *img) {
630  int ret=0;
631 
632  if(IMG_TEST) printf("\nimgReadAnalyzeFirstFrame(%s, *img)\n", fname);
633  /* Check the input */
634  if(img==NULL) return STATUS_FAULT;
635  if(img->status!=IMG_STATUS_INITIALIZED) return STATUS_FAULT;
637  if(fname==NULL) return STATUS_FAULT;
638 
639  /* Read header information from file */
640  ret=imgReadAnalyzeHeader(fname, img); if(ret) return(ret);
641  if(IMG_TEST>3) imgInfo(img);
642 
643  /* Allocate memory for one frame */
644  img->dimt=1;
645  ret=imgAllocate(img, img->dimz, img->dimy, img->dimx, img->dimt);
646  if(ret) return STATUS_NOMEMORY;
647 
648  /* Read the first frame */
649  ret=imgReadAnalyzeFrame(fname, 1, img, 0); if(ret) return(ret);
650 
651  /* All went well */
652  imgSetStatus(img, STATUS_OK);
653  return STATUS_OK;
654 }
655 /*****************************************************************************/
656 
657 /*****************************************************************************/
674 int imgReadAnalyzeFrame(const char *fname, int frame_to_read, IMG *img, int frame_index) {
675  FILE *fp;
676  int ret, pi, xi, yi;
677  float *fdata=NULL, *fptr;
678  ANALYZE_DSR dsr;
679  char datfile[FILENAME_MAX], hdrfile[FILENAME_MAX], siffile[FILENAME_MAX];
680  SIF sif;
681 
682 
683  if(IMG_TEST) printf("\nimgReadAnalyzeFrame(%s, %d, *img, %d)\n",
684  fname, frame_to_read, frame_index);
685 
686  /* Check the input */
687  if(img==NULL) return STATUS_FAULT;
688  if(img->status!=IMG_STATUS_OCCUPIED) return STATUS_FAULT;
689  if(fname==NULL) return STATUS_FAULT;
690  if(frame_index<0 || frame_index>img->dimt-1) return STATUS_FAULT;
691  if(frame_to_read<1) return STATUS_FAULT;
693 
694  /* Determine the names of hdr, data and sif files */
695  ret=anaDatabaseExists(fname, hdrfile, datfile, siffile);
696  if(ret==0) return STATUS_NOFILE;
697 
698  /* Read Analyze header file */
699  ret=anaReadHeader(hdrfile, &dsr);
700  if(ret!=0) {
701  if(ret==1) return STATUS_FAULT;
702  else if(ret==2) return STATUS_NOHEADERFILE;
703  else return STATUS_UNSUPPORTED;
704  return(STATUS_FAULT);
705  }
706 
707  /* Open image datafile */
708  if(IMG_TEST>2) fprintf(stdout, "reading image data %s\n", datfile);
709  if((fp=fopen(datfile, "rb")) == NULL) return STATUS_NOIMGDATA;
710 
711  /* Allocate memory for one image frame */
712  fdata=malloc(img->dimx*img->dimy*img->dimz*sizeof(float));
713  if(fdata==NULL) {fclose(fp); return STATUS_NOMEMORY;}
714 
715  /* Read the required image frame */
716  fptr=fdata;
717  ret=anaReadImagedata(fp, &dsr, frame_to_read, fptr);
718  fclose(fp);
719  if(ret==3) {free(fdata); return STATUS_NOMATRIX;} /* no more frames */
720  if(ret!=0) {free(fdata); return STATUS_UNSUPPORTED;}
721 
722  /* Copy pixel values to IMG */
723  fptr=fdata;
724  if(anaFlipping()==0) { /* no flipping in z-direction */
725  for(pi=0; pi<img->dimz; pi++)
726  for(yi=img->dimy-1; yi>=0; yi--)
727  for(xi=img->dimx-1; xi>=0; xi--)
728  img->m[pi][yi][xi][frame_index]=*fptr++;
729  } else {
730  for(pi=img->dimz-1; pi>=0; pi--)
731  for(yi=img->dimy-1; yi>=0; yi--)
732  for(xi=img->dimx-1; xi>=0; xi--)
733  img->m[pi][yi][xi][frame_index]=*fptr++;
734  }
735  free(fdata);
736 
737  imgSetStatus(img, STATUS_OK); /* If the rest is failed, no problem */
738 
739  /*
740  * Try to read frame time information from SIF file
741  */
742  sifInit(&sif);
743  if(sifRead(siffile, &sif)!=0) return STATUS_OK;
744  /* Frame information */
745  if(sif.frameNr>=frame_to_read) {
746  img->start[frame_index]=sif.x1[frame_to_read-1];
747  img->end[frame_index]=sif.x2[frame_to_read-1];
748  img->mid[frame_index]=0.5*(img->start[frame_index]+img->end[frame_index]);
749  img->prompts[frame_index]=sif.prompts[frame_to_read-1];
750  img->randoms[frame_index]=sif.randoms[frame_to_read-1];
751  }
752  sifEmpty(&sif);
753 
754  return STATUS_OK;
755 }
756 /*****************************************************************************/
757 
758 /*****************************************************************************/
781 int imgWriteAnalyzeFrame(const char *dbname, int frame_to_write, IMG *img, int frame_index, float fmin, float fmax) {
782  IMG test_img;
783  int ret=0, pxlNr, zi, xi, yi, little;
784  FILE *fp;
785  short int *sdata=NULL, *sptr;
786  char datfile[FILENAME_MAX], hdrfile[FILENAME_MAX], siffile[FILENAME_MAX];
787  ANALYZE_DSR dsr;
788  float scale_factor=1.0;
789 
790 
791  if(IMG_TEST) printf("\nimgWriteAnalyzeFrame(%s, %d, *img, %d, %g, %g)\n",
792  dbname, frame_to_write, frame_index, fmin, fmax);
793 
794  /*
795  * Check the input
796  */
797  if(dbname==NULL) return STATUS_FAULT;
798  if(img==NULL) return STATUS_FAULT;
799  if(img->status!=IMG_STATUS_OCCUPIED) return STATUS_FAULT;
800  if(frame_to_write<0) return STATUS_FAULT;
801  if(frame_index<0 || frame_index>=img->dimt) return STATUS_FAULT;
802  if(img->_fileFormat!=IMG_ANA_L && img->_fileFormat!=IMG_ANA) return STATUS_FAULT;
803 
804  /*
805  * If database does not exist, then create it with new header,
806  * and if it does exist, then read and check header information.
807  * Create or edit header to contain correct frame nr.
808  * Determine the global scaling factor.
809  */
810  imgInit(&test_img);
811  if(anaDatabaseExists(dbname, hdrfile, datfile, siffile)==0) { /* not existing */
812 
813  /* Create database filenames */
814  sprintf(hdrfile, "%s.hdr", dbname);
815  sprintf(datfile, "%s.img", dbname);
816  sprintf(siffile, "%s.sif", dbname);
817 
818  /* Set main header */
819  imgSetAnalyzeHeader(img, dbname, &dsr, fmin, fmax);
820  if(frame_to_write==0) frame_to_write=1;
821  dsr.dime.dim[4]=frame_to_write;
822  scale_factor=dsr.dime.funused1;
823  if(fabs(scale_factor)>1.0E-20) scale_factor=1.0/scale_factor;
824 
825  /* Write Analyze header */
826  ret=anaWriteHeader(hdrfile, &dsr); if(ret && IMG_TEST) printf("anaWriteHeader() := %d\n", ret);
827  if(ret) return STATUS_CANTWRITEHEADERFILE;
828 
829  /* Remove datafile if necessary */
830  if(access(datfile, 0) != -1) remove(datfile);
831 
832  } else { /* database does exist */
833 
834  /* Read header information for checking */
835  ret=imgReadAnalyzeHeader(dbname, &test_img); if(ret!=0) return ret;
836  /* Check that file format is the same */
837  if(img->_fileFormat!=test_img._fileFormat || img->type!=test_img.type)
838  return STATUS_WRONGFILETYPE;
839  /* Check that matrix sizes are the same */
840  if(img->dimz!=test_img.dimz || img->dimx!=test_img.dimx ||
841  img->dimy!=test_img.dimy)
842  return STATUS_VARMATSIZE;
843  imgEmpty(&test_img);
844 
845  /* Read the header, set new frame number, and write it back */
846  /* Get also the scale factor */
847  if((ret=anaReadHeader(hdrfile, &dsr))!=0) return STATUS_NOMAINHEADER;
848  scale_factor=1.0/dsr.dime.funused1;
849  if(frame_to_write==0) frame_to_write=dsr.dime.dim[4]+1;
850  if(dsr.dime.dim[4]<frame_to_write) {
851  if(dsr.dime.dim[4]+1<frame_to_write) return STATUS_MISSINGMATRIX;
852  dsr.dime.dim[4]=frame_to_write;
853  }
854  if((ret=anaWriteHeader(hdrfile, &dsr))!=0) return STATUS_NOWRITEPERM;
855  }
856  if(IMG_TEST>2) {
857  printf("frame_to_write := %d\n", frame_to_write);
858  printf("hdrfile := %s\n", hdrfile);
859  printf("datfile := %s\n", datfile);
860  printf("siffile := %s\n", siffile);
861  }
862 
863  /* Allocate memory for matrix short int data (one plane) */
864  pxlNr=img->dimx*img->dimy;
865  sdata=(short int*)malloc(pxlNr*sizeof(short int));
866  if(sdata==NULL) return STATUS_NOMEMORY;
867 
868  /* Open datafile, not removing possible old contents */
869  if(frame_to_write==1) fp=fopen(datfile, "wb"); else fp=fopen(datfile, "r+b");
870  if(fp==NULL) {free(sdata); return STATUS_CANTWRITEIMGFILE;}
871  /* Move file pointer to the place of current frame */
872  if(fseek(fp, (frame_to_write-1)*pxlNr*img->dimz*sizeof(short int), SEEK_SET)!=0) {
873  free(sdata); fclose(fp); return STATUS_MISSINGMATRIX;}
874  little=little_endian();
875  /* Copy, scale and write data plane-by-plane */
876  if(anaFlipping()==0) {
877  for(zi=0; zi<img->dimz; zi++) {
878  sptr=sdata; /*printf("plane := %d\n scale_factor := %g\n", zi+1, scale_factor);*/
879  for(yi=img->dimy-1; yi>=0; yi--) for(xi=img->dimx-1; xi>=0; xi--) {
880  *sptr=temp_roundf(scale_factor*img->m[zi][yi][xi][frame_index]); sptr++;
881  }
882  /* Change byte order if necessary */
883  sptr=sdata; if(little!=dsr.little) swabip(sptr, pxlNr*sizeof(short int));
884  /* Write image data */
885  sptr=sdata;
886  if(fwrite(sptr, sizeof(short int), pxlNr, fp) != pxlNr) {
887  free(sdata); fclose(fp); return STATUS_CANTWRITEIMGFILE;
888  }
889  }
890  } else {
891  for(zi=img->dimz-1; zi>=0; zi--) {
892  sptr=sdata;
893  for(yi=img->dimy-1; yi>=0; yi--) for(xi=img->dimx-1; xi>=0; xi--) {
894  *sptr=temp_roundf(scale_factor*img->m[zi][yi][xi][frame_index]); sptr++;
895  }
896  /* Change byte order if necessary */
897  sptr=sdata; if(little!=dsr.little) swabip(sptr, pxlNr*sizeof(short int));
898  /* Write image data */
899  sptr=sdata;
900  if(fwrite(sptr, sizeof(short int), pxlNr, fp) != pxlNr) {
901  free(sdata); fclose(fp); return STATUS_CANTWRITEIMGFILE;
902  }
903  }
904  }
905  free(sdata);
906  fclose(fp);
907 
908  return STATUS_OK;
909 }
910 /*****************************************************************************/
911 
912 /*****************************************************************************/
int anaExists(const char *dbname)
Definition: analyze.c:75
time_t scantime
Definition: sif.h:38
int imgAllocate(IMG *image, int planes, int rows, int columns, int frames)
Definition: img.c:285
int anaWriteHeader(char *filename, ANALYZE_DSR *h)
Definition: analyze.c:209
void imgEmpty(IMG *image)
Definition: img.c:216
unsigned short int dimz
Definition: img.h:265
void sifInit(SIF *data)
Definition: sif.c:61
Definition: img.h:156
char patient_id[10]
Definition: analyze.h:85
char db_name[18]
Definition: analyze.h:46
ANALYZE_HEADER_HISTORY hist
Definition: analyze.h:102
float sizey
Definition: img.h:210
short int datatype
Definition: analyze.h:62
float * end
Definition: img.h:292
int anaFlipping()
Definition: analyze.c:545
unsigned short int dimx
Definition: img.h:261
int anaReadHeader(char *filename, ANALYZE_DSR *h)
Definition: analyze.c:102
int imgReadAnalyzeFrame(const char *fname, int frame_to_read, IMG *img, int frame_index)
Definition: img_ana.c:674
char patientName[32]
Definition: img.h:176
int imgWriteAnalyze(const char *dbname, IMG *img)
Definition: img_ana.c:253
float * mid
Definition: img.h:294
int frameNr
Definition: sif.h:40
char status
Definition: img.h:164
float **** m
Definition: img.h:274
double * prompts
Definition: sif.h:54
#define IMG_STATUS_OCCUPIED
Definition: img.h:73
float * randoms
Definition: img.h:308
short int dim[8]
Definition: analyze.h:54
int imgSetAnalyzeHeader(IMG *img, const char *dbname, ANALYZE_DSR *dsr, float fmin, float fmax)
Definition: img_ana.c:542
#define IMG_ANA_L
Definition: img.h:90
int anaPrintHeader(ANALYZE_DSR *h, FILE *fp)
Definition: analyze.c:307
float sizez
Definition: img.h:212
Definition: img.h:118
char data_type[10]
Definition: analyze.h:45
int little
Definition: analyze.h:104
double * randoms
Definition: sif.h:56
int sifRead(char *filename, SIF *data)
Definition: sifio.c:64
int imgReadAnalyze(const char *dbname, IMG *img)
Definition: img_ana.c:83
short int bitpix
Definition: analyze.h:63
#define IMG_STATUS_INITIALIZED
Definition: img.h:72
Definition: sif.h:36
double * x2
Definition: sif.h:52
char studynr[11]
Definition: sif.h:46
#define IMG_ANA
Definition: img.h:89
float isotopeHalflife
Definition: img.h:182
int anaDatabaseExists(const char *dbname, char *hdrfile, char *imgfile, char *siffile)
Definition: analyze.c:619
void imgInit(IMG *image)
Definition: img.c:163
int imgReadAnalyzeFirstFrame(const char *fname, IMG *img)
Definition: img_ana.c:629
int anaReadImagedata(FILE *fp, ANALYZE_DSR *h, int frame, float *data)
Definition: analyze.c:381
#define IMG_TYPE_IMAGE
Definition: img.h:80
#define ANALYZE_DT_SIGNED_SHORT
Definition: analyze.h:33
int _fileFormat
Definition: img.h:229
float sizex
Definition: img.h:208
char studyNr[MAX_STUDYNR_LEN+1]
Definition: img.h:174
int imgReadAnalyzeHeader(const char *dbname, IMG *img)
Definition: img_ana.c:414
char decayCorrected
Definition: img.h:184
double * x1
Definition: sif.h:50
int IMG_TEST
Definition: img.h:128
int imgWriteAnalyzeFrame(const char *dbname, int frame_to_write, IMG *img, int frame_index, float fmin, float fmax)
Definition: img_ana.c:781
time_t scanStart
Definition: img.h:186
char isotope_name[8]
Definition: sif.h:48
ANALYZE_HEADER_IMGDIM dime
Definition: analyze.h:101
int * planeNumber
Definition: img.h:284
float * start
Definition: img.h:290
unsigned short int dimy
Definition: img.h:263
unsigned short int dimt
Definition: img.h:259
void imgInfo(IMG *image)
Definition: img.c:414
ANALYZE_HEADER_KEY hk
Definition: analyze.h:100
float * prompts
Definition: img.h:306
void sifEmpty(SIF *data)
Definition: sif.c:74
int imgGetAnalyzeHeader(IMG *img, ANALYZE_DSR *h)
Definition: img_ana.c:481
char type
Definition: img.h:198
void imgSetStatus(IMG *img, int status_index)
Definition: img.c:399
int imgMinMax(IMG *img, float *minvalue, float *maxvalue)
Definition: imgmax.c:115