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/*****************************************************************************/
77int 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/*****************************************************************************/
374int 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;
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/*****************************************************************************/
568int 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/*****************************************************************************/
886int 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;
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;
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/*****************************************************************************/
1289int 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;
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/*****************************************************************************/
1422int 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;
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/*****************************************************************************/
1464int 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;
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);
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/*****************************************************************************/
1656int 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);
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/*****************************************************************************/
1790void 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
int ecat63WriteImage(FILE *fp, int matnum, ECAT63_imageheader *h, void *data)
Definition ecat63w.c:429
FILE * ecat63Create(const char *fname, ECAT63_mainheader *h)
Definition ecat63w.c:386
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
int matstat
Definition ecat63.h:59
int endblk
Definition ecat63.h:58
int matnum
Definition ecat63.h:56
int strtblk
Definition ecat63.h:57
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
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 plane_eff_corr_fctr
Definition ecat63.h:117
short int dimension_1
Definition ecat63.h:107
float ecat_calibration_fctr
Definition ecat63.h:121
float loss_corr_fctr
Definition ecat63.h:118
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
float loss_correction_fctr
Definition ecat63.h:142
short int dimension_2
Definition ecat63.h:129
short int scan_max
Definition ecat63.h:136
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