My Project
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/*****************************************************************************/
83int 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 */
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);
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/*****************************************************************************/
253int 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 */
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 */
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 }
399 return 0;
400}
401/*****************************************************************************/
402
403/*****************************************************************************/
414int 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;
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]) {
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
525 return STATUS_OK;
526}
527/*****************************************************************************/
528
529/*****************************************************************************/
542int 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
616 return STATUS_OK;
617}
618/*****************************************************************************/
619
620/*****************************************************************************/
629int 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;
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 */
653 return STATUS_OK;
654}
655/*****************************************************************************/
656
657/*****************************************************************************/
674int 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/*****************************************************************************/
781int 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)
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 anaWriteHeader(char *filename, ANALYZE_DSR *h)
Definition: analyze.c:210
int anaFlipping()
Definition: analyze.c:546
int anaPrintHeader(ANALYZE_DSR *h, FILE *fp)
Definition: analyze.c:308
int anaDatabaseExists(const char *dbname, char *hdrfile, char *imgfile, char *siffile)
Definition: analyze.c:620
int anaExists(const char *dbname)
Definition: analyze.c:76
int anaReadHeader(char *filename, ANALYZE_DSR *h)
Definition: analyze.c:103
int anaReadImagedata(FILE *fp, ANALYZE_DSR *h, int frame, float *data)
Definition: analyze.c:382
#define ANALYZE_DT_SIGNED_SHORT
Definition: analyze.h:33
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
@ STATUS_NOHEADERFILE
Definition: img.h:122
@ STATUS_NOSIFDATA
Definition: img.h:123
@ STATUS_NOIMGDATA
Definition: img.h:123
@ STATUS_WRONGSIFDATA
Definition: img.h:123
@ STATUS_INVALIDHEADER
Definition: img.h:122
@ STATUS_OK
Definition: img.h:118
@ STATUS_WRONGFILETYPE
Definition: img.h:124
@ STATUS_NOWRITEPERM
Definition: img.h:119
@ STATUS_NOFILE
Definition: img.h:118
@ STATUS_CANTWRITEIMGFILE
Definition: img.h:124
@ STATUS_NOMATRIX
Definition: img.h:121
@ STATUS_CANTWRITEHEADERFILE
Definition: img.h:124
@ STATUS_MISSINGMATRIX
Definition: img.h:119
@ 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_ANA_L
Definition: img.h:90
#define IMG_STATUS_INITIALIZED
Definition: img.h:72
#define IMG_ANA
Definition: img.h:89
#define IMG_TYPE_IMAGE
Definition: img.h:80
int imgWriteAnalyze(const char *dbname, IMG *img)
Definition: img_ana.c:253
int imgGetAnalyzeHeader(IMG *img, ANALYZE_DSR *h)
Definition: img_ana.c:481
int imgReadAnalyzeFirstFrame(const char *fname, IMG *img)
Definition: img_ana.c:629
int imgReadAnalyzeFrame(const char *fname, int frame_to_read, IMG *img, int frame_index)
Definition: img_ana.c:674
int imgReadAnalyzeHeader(const char *dbname, IMG *img)
Definition: img_ana.c:414
int imgSetAnalyzeHeader(IMG *img, const char *dbname, ANALYZE_DSR *dsr, float fmin, float fmax)
Definition: img_ana.c:542
int imgWriteAnalyzeFrame(const char *dbname, int frame_to_write, IMG *img, int frame_index, float fmin, float fmax)
Definition: img_ana.c:781
int imgReadAnalyze(const char *dbname, IMG *img)
Definition: img_ana.c:83
int imgMinMax(IMG *img, float *minvalue, float *maxvalue)
Definition: imgmax.c:115
void sifInit(SIF *data)
Definition: sif.c:61
void sifEmpty(SIF *data)
Definition: sif.c:74
int sifRead(char *filename, SIF *data)
Definition: sifio.c:64
ANALYZE_HEADER_HISTORY hist
Definition: analyze.h:102
int little
Definition: analyze.h:104
ANALYZE_HEADER_KEY hk
Definition: analyze.h:100
ANALYZE_HEADER_IMGDIM dime
Definition: analyze.h:101
char patient_id[10]
Definition: analyze.h:85
short int datatype
Definition: analyze.h:62
short int dim[8]
Definition: analyze.h:54
short int bitpix
Definition: analyze.h:63
char db_name[18]
Definition: analyze.h:46
char data_type[10]
Definition: analyze.h:45
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
float **** m
Definition: img.h:274
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
unsigned short int dimt
Definition: img.h:259
int * planeNumber
Definition: img.h:284
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 isotopeHalflife
Definition: img.h:182
char studyNr[MAX_STUDYNR_LEN+1]
Definition: img.h:174
float * randoms
Definition: img.h:308
float * mid
Definition: img.h:294
float sizez
Definition: img.h:212
Definition: sif.h:36
double * x1
Definition: sif.h:50
double * prompts
Definition: sif.h:54
char studynr[11]
Definition: sif.h:46
int frameNr
Definition: sif.h:40
double * x2
Definition: sif.h:52
time_t scantime
Definition: sif.h:38
char isotope_name[8]
Definition: sif.h:48
double * randoms
Definition: sif.h:56