00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044 #include <stdio.h>
00045 #include <stdlib.h>
00046 #include <unistd.h>
00047 #include <math.h>
00048 #include <string.h>
00049 #include <time.h>
00050
00051 #include "petc99.h"
00052 #include "swap.h"
00053 #include "halflife.h"
00054
00055 #include "include/img.h"
00056 #include "include/ecat63.h"
00057 #include "include/ecat7.h"
00058 #include "include/imgmax.h"
00059 #include "include/imgdecay.h"
00060 #include "include/sif.h"
00061 #include "include/imgfile.h"
00062
00063
00064
00077 int ecat63ReadAllToImg(const char *fname, IMG *img) {
00078 int i, j, m, ret, blkNr=0, dim_x=0, dim_y=0, pxlNr=0;
00079 int frameNr, planeNr, del_nr=0;
00080 int frame, plane, prev_frame, prev_plane, seqplane;
00081 FILE *fp;
00082 ECAT63_mainheader main_header;
00083 MATRIXLIST mlist;
00084 MatDir tmpMatdir;
00085 Matval matval;
00086 ECAT63_imageheader image_header;
00087 ECAT63_scanheader scan_header;
00088 ECAT63_attnheader attn_header;
00089 ECAT63_normheader norm_header;
00090 float scale=1.0;
00091 short int *sptr;
00092 char *mdata=NULL, *mptr;
00093 int *iptr;
00094 struct tm scanStart;
00095
00096
00097 if(IMG_TEST) printf("ecat63ReadAllToImg(%s, *img)\n", fname);
00098
00099 if(fname==NULL || img==NULL || img->status!=IMG_STATUS_INITIALIZED) {
00100 if(img) imgSetStatus(img, STATUS_FAULT);
00101 return 1;
00102 }
00103
00104
00105 if((fp=fopen(fname, "rb")) == NULL) {
00106 imgSetStatus(img, STATUS_NOFILE);
00107 return 3;
00108 }
00109
00110
00111 if((ret=ecat63ReadMainheader(fp, &main_header))) {
00112 imgSetStatus(img, STATUS_NOMAINHEADER);
00113 return 4;
00114 }
00115 if(IMG_TEST>5) ecat63PrintMainheader(&main_header, stdout);
00116
00117
00118 ecat63InitMatlist(&mlist);
00119 ret=ecat63ReadMatlist(fp, &mlist);
00120 if(ret) {
00121 imgSetStatus(img, STATUS_NOMATLIST);
00122 return 5;
00123 }
00124 if(mlist.matrixNr<=0) {
00125 strcpy(ecat63errmsg, "matrix list is empty");
00126 return 5;
00127 }
00128
00129 for(i=0; i<mlist.matrixNr-1; i++) for(j=i+1; j<mlist.matrixNr; j++) {
00130 if(mlist.matdir[i].matnum>mlist.matdir[j].matnum) {
00131 tmpMatdir=mlist.matdir[i];
00132 mlist.matdir[i]=mlist.matdir[j]; mlist.matdir[j]=tmpMatdir;
00133 }
00134 }
00135 if(IMG_TEST>6) {
00136 printf("matrix list after sorting:\n");
00137 ecat63PrintMatlist(&mlist);
00138 }
00139
00140
00141 if(main_header.num_frames>0) {
00142 del_nr=ecat63DeleteLateFrames(&mlist, main_header.num_frames);
00143 if(IMG_TEST && del_nr>0)
00144 printf(" %d entries in matrix list will not be used.\n", del_nr);
00145 }
00146
00147
00148
00149 prev_plane=plane=-1; prev_frame=frame=-1; frameNr=planeNr=0; ret=0;
00150 for(m=0; m<mlist.matrixNr; m++) if(mlist.matdir[m].matstat==1) {
00151
00152 mat_numdoc(mlist.matdir[m].matnum, &matval);
00153 plane=matval.plane;
00154 if(main_header.num_frames>=main_header.num_gates) frame=matval.frame;
00155 else frame=matval.gate;
00156 if(IMG_TEST>2) printf("frame=%d plane=%d\n", frame, plane);
00157
00158 if(plane!=prev_plane) {
00159 frameNr=1; planeNr++;
00160 } else {
00161 frameNr++;
00162
00163 if(prev_frame>0 && frame!=prev_frame+1) {ret=1; break;}
00164 }
00165 prev_plane=plane; prev_frame=frame;
00166
00167 if(m==0) {
00168 blkNr=mlist.matdir[m].endblk-mlist.matdir[m].strtblk;
00169 } else if(blkNr!=mlist.matdir[m].endblk-mlist.matdir[m].strtblk) {
00170 ret=2; break;
00171 }
00172 }
00173 if(IMG_TEST) printf("frameNr=%d planeNr=%d\n", frameNr, planeNr);
00174 if(ret==1 || (frameNr*planeNr != mlist.matrixNr-del_nr)) {
00175 strcpy(ecat63errmsg, "missing matrix");
00176 ecat63EmptyMatlist(&mlist); fclose(fp);
00177 return(6);
00178 }
00179 if(ret==2) {
00180 strcpy(ecat63errmsg, "matrix sizes are different");
00181 ecat63EmptyMatlist(&mlist); fclose(fp); return(7);
00182 }
00183
00184 m=0; if(main_header.file_type==IMAGE_DATA) {
00185 ret=ecat63ReadImageheader(fp, mlist.matdir[m].strtblk, &image_header);
00186 dim_x=image_header.dimension_1; dim_y=image_header.dimension_2;
00187 } else if(main_header.file_type==RAW_DATA) {
00188 ret=ecat63ReadScanheader(fp, mlist.matdir[m].strtblk, &scan_header);
00189 dim_x=scan_header.dimension_1; dim_y=scan_header.dimension_2;
00190 } else if(main_header.file_type==ATTN_DATA) {
00191 ret=ecat63ReadAttnheader(fp, mlist.matdir[m].strtblk, &attn_header);
00192 dim_x=attn_header.dimension_1; dim_y=attn_header.dimension_2;
00193 } else if(main_header.file_type==NORM_DATA) {
00194 ret=ecat63ReadNormheader(fp, mlist.matdir[m].strtblk, &norm_header);
00195 dim_x=norm_header.dimension_1; dim_y=norm_header.dimension_2;
00196 }
00197 pxlNr=dim_x*dim_y;
00198 if(ret) {
00199 sprintf(ecat63errmsg, "cannot read matrix %u subheader", mlist.matdir[m].matnum);
00200 ecat63EmptyMatlist(&mlist); fclose(fp); return(8);
00201 }
00202
00203
00204 if(IMG_TEST>1) printf("allocating memory for %d blocks\n", blkNr);
00205 mdata=(char*)malloc(blkNr*MatBLKSIZE); if(mdata==NULL) {
00206 strcpy(ecat63errmsg, "out of memory");
00207 fclose(fp); ecat63EmptyMatlist(&mlist); return(8);
00208 }
00209
00210 ret=imgAllocate(img, planeNr, dim_y, dim_x, frameNr);
00211 if(ret) {
00212 sprintf(ecat63errmsg, "out of memory (%d)", ret);
00213 fclose(fp); ecat63EmptyMatlist(&mlist); free(mdata); return(9);
00214 }
00215
00216
00217 img->scanner=main_header.system_type;
00218 img->unit=main_header.calibration_units;
00219 strncpy(img->radiopharmaceutical, main_header.radiopharmaceutical, 32);
00220 img->isotopeHalflife=main_header.isotope_halflife;
00221 memset(&scanStart, 0, sizeof(struct tm));
00222 scanStart.tm_year=main_header.scan_start_year-1900;
00223 scanStart.tm_mon=main_header.scan_start_month-1;
00224 scanStart.tm_mday=main_header.scan_start_day; scanStart.tm_yday=0;
00225 scanStart.tm_hour=main_header.scan_start_hour;
00226 scanStart.tm_min=main_header.scan_start_minute;
00227 scanStart.tm_sec=main_header.scan_start_second;
00228 scanStart.tm_isdst=-1;
00229 img->scanStart=mktime(&scanStart); if(img->scanStart==-1) img->scanStart=0;
00230 img->axialFOV=10.*main_header.axial_fov;
00231 img->transaxialFOV=10.*main_header.transaxial_fov;
00232 strncpy(img->studyNr, main_header.study_name, MAX_STUDYNR_LEN);
00233 strncpy(img->patientName, main_header.patient_name, 31);
00234 strncpy(img->patientID, main_header.patient_id, 15);
00235 img->sizez=10.*main_header.plane_separation;
00236 if(main_header.file_type==IMAGE_DATA) {
00237 img->type=IMG_TYPE_IMAGE;
00238 if(img->unit<1) img->unit=image_header.quant_units;
00239 img->zoom=image_header.recon_scale;
00240 if(image_header.decay_corr_fctr>1.0) img->decayCorrected=1;
00241 img->sizex=img->sizey=10.*image_header.pixel_size;
00242 if(img->sizez<10.*image_header.slice_width)
00243 img->sizez=10.*image_header.slice_width;
00244 } else if(main_header.file_type==RAW_DATA) {
00245 img->type=IMG_TYPE_RAW;
00246 img->decayCorrected=0;
00247 } else if(main_header.file_type==ATTN_DATA) {
00248 img->type=IMG_TYPE_RAW;
00249 img->decayCorrected=0;
00250 } else if(main_header.file_type==NORM_DATA) {
00251 img->type=IMG_TYPE_RAW;
00252 img->decayCorrected=0;
00253 }
00254 strncpy(img->studyDescription, main_header.study_description, 32);
00255 strncpy(img->userProcessCode, main_header.user_process_code, 10);
00256 img->userProcessCode[10]=(char)0;
00257
00258 if(!img->studyNr[0] && studynr_validity_check(img->userProcessCode))
00259 strcpy(img->studyNr, img->userProcessCode);
00260
00261
00262 img->_fileFormat=IMG_E63;
00263
00264
00265 prev_plane=plane=-1; seqplane=-1;
00266 for(m=0; m<mlist.matrixNr; m++) if(mlist.matdir[m].matstat==1) {
00267
00268 mat_numdoc(mlist.matdir[m].matnum, &matval);
00269 plane=matval.plane;
00270
00271 if(plane!=prev_plane) {seqplane++; frame=1;} else frame++;
00272 prev_plane=plane;
00273 img->planeNumber[seqplane]=plane;
00274
00275 if(main_header.file_type==IMAGE_DATA) {
00276 ret=ecat63ReadImageheader(fp, mlist.matdir[m].strtblk, &image_header);
00277 if(dim_x!=image_header.dimension_1 || dim_y!=image_header.dimension_2) ret=1;
00278 scale=image_header.quant_scale;
00279 if(image_header.ecat_calibration_fctr>0.0
00280 && fabs(main_header.calibration_factor/image_header.ecat_calibration_fctr-1.0)>0.0001)
00281 scale*=image_header.ecat_calibration_fctr;
00282 if(img->unit==0 && image_header.quant_units>0) img->unit=image_header.quant_units;
00283 img->_dataType=image_header.data_type;
00284 img->start[frame-1]=image_header.frame_start_time/1000.;
00285 img->end[frame-1]=img->start[frame-1]+image_header.frame_duration/1000.;
00286 img->mid[frame-1]=0.5*(img->start[frame-1]+img->end[frame-1]);
00287 if(image_header.decay_corr_fctr>1.0)
00288 img->decayCorrFactor[frame-1]=image_header.decay_corr_fctr;
00289 } else if(main_header.file_type==RAW_DATA) {
00290 ret=ecat63ReadScanheader(fp, mlist.matdir[m].strtblk, &scan_header);
00291 if(dim_x!=scan_header.dimension_1 || dim_y!=scan_header.dimension_2) ret=1;
00292 scale=scan_header.scale_factor;
00293 if(scan_header.loss_correction_fctr>0.0) scale*=scan_header.loss_correction_fctr;
00294 img->_dataType=scan_header.data_type;
00295 img->start[frame-1]=scan_header.frame_start_time/1000.;
00296 img->end[frame-1]=img->start[frame-1]+scan_header.frame_duration/1000.;
00297 img->mid[frame-1]=0.5*(img->start[frame-1]+img->end[frame-1]);
00298 img->sampleDistance=10.0*scan_header.sample_distance;
00299 img->prompts[frame-1]=(float)scan_header.prompts;
00300 img->randoms[frame-1]=(float)scan_header.delayed;
00301 } else if(main_header.file_type==ATTN_DATA) {
00302 ret=ecat63ReadAttnheader(fp, mlist.matdir[m].strtblk, &attn_header);
00303 if(dim_x!=attn_header.dimension_1 || dim_y!=attn_header.dimension_2) ret=1;
00304 scale=attn_header.scale_factor;
00305 img->_dataType=attn_header.data_type;
00306 } else if(main_header.file_type==NORM_DATA) {
00307 ret=ecat63ReadNormheader(fp, mlist.matdir[m].strtblk, &norm_header);
00308 if(dim_x!=norm_header.dimension_1 || dim_y!=norm_header.dimension_2) ret=1;
00309 scale=norm_header.scale_factor;
00310 img->_dataType=norm_header.data_type;
00311 } else img->_dataType=-1;
00312 if(ret) {
00313 sprintf(ecat63errmsg, "cannot read matrix %u subheader", mlist.matdir[m].matnum);
00314 ecat63EmptyMatlist(&mlist); free(mdata); fclose(fp); return(10);
00315 }
00316 if(main_header.calibration_factor>0.0) scale*=main_header.calibration_factor;
00317 if(IMG_TEST>2) printf("scale=%e datatype=%d\n", scale, img->_dataType);
00318
00319 ret=ecat63ReadMatdata(fp, mlist.matdir[m].strtblk+1,
00320 mlist.matdir[m].endblk-mlist.matdir[m].strtblk,
00321 mdata, img->_dataType);
00322 if(ret) {
00323 strcpy(ecat63errmsg, "cannot read matrix data");
00324 ecat63EmptyMatlist(&mlist); free(mdata); fclose(fp); return(11);
00325 }
00326 if(img->_dataType==BYTE_TYPE) {
00327 for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr++) {
00328 img->m[seqplane][i][j][frame-1]=scale*(float)(*mptr);
00329 }
00330 } else if(img->_dataType==VAX_I2 || img->_dataType==SUN_I2) {
00331 for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr+=2) {
00332 sptr=(short int*)mptr;
00333 img->m[seqplane][i][j][frame-1]=scale*(float)(*sptr);
00334 }
00335 } else if(img->_dataType==VAX_I4 || img->_dataType==SUN_I4) {
00336 for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr+=4) {
00337 iptr=(int*)mptr;
00338 img->m[seqplane][i][j][frame-1]=1.0;
00339 }
00340 } else if(img->_dataType==VAX_R4 || img->_dataType==IEEE_R4) {
00341 for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr+=4) {
00342 memcpy(&img->m[seqplane][i][j][frame-1], mptr, 4);
00343 img->m[seqplane][i][j][frame-1]*=scale;
00344 }
00345 }
00346 }
00347
00348
00349 imgUnitFromEcat(img, img->unit);
00350
00351
00352 free(mdata);
00353 ecat63EmptyMatlist(&mlist);
00354 fclose(fp);
00355
00356
00357 if(img->_dataType==VAX_I4 || img->_dataType==VAX_R4) img->_dataType=VAX_I2;
00358
00359 return(0);
00360 }
00361
00362
00363
00374 int ecat63WriteAllImg(const char *fname, IMG *img) {
00375 int frame, plane, n, i, j, m, ret=0;
00376 float f, fmax, fmin, g, scale;
00377 short int *sdata, *sptr, smin, smax;
00378 FILE *fp;
00379 ECAT63_mainheader main_header;
00380 ECAT63_imageheader image_header;
00381 ECAT63_scanheader scan_header;
00382 struct tm *scanStart;
00383
00384
00385 if(IMG_TEST) printf("ecat63WriteAllImg(%s, *img)\n", fname);
00386
00387 if(fname==NULL) {strcpy(ecat63errmsg, "invalid ECAT filename"); return(1);}
00388 if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) {
00389 strcpy(ecat63errmsg, "image data is empty"); return(2);}
00390 if(img->_dataType<1) img->_dataType=VAX_I2;
00391
00392
00393 memset(&main_header, 0, sizeof(ECAT63_mainheader));
00394 memset(&image_header,0, sizeof(ECAT63_imageheader));
00395 memset(&scan_header, 0, sizeof(ECAT63_scanheader));
00396
00397
00398 main_header.sw_version=2;
00399 main_header.num_planes=img->dimz;
00400 main_header.num_frames=img->dimt;
00401 main_header.num_gates=1;
00402 main_header.num_bed_pos=1;
00403 if(img->type==IMG_TYPE_IMAGE) main_header.file_type=IMAGE_DATA;
00404 else main_header.file_type=RAW_DATA;
00405 main_header.data_type=img->_dataType;
00406 if(img->scanner>0) main_header.system_type=img->scanner;
00407 else main_header.system_type=ECAT63_SYSTEM_TYPE_DEFAULT;
00408 main_header.calibration_factor=1.0;
00409 main_header.axial_fov=img->axialFOV/10.0;
00410 main_header.transaxial_fov=img->transaxialFOV/10.0;
00411 main_header.plane_separation=img->sizez/10.0;
00412 main_header.calibration_units=imgUnitToEcat6(img);
00413 strncpy(main_header.radiopharmaceutical, img->radiopharmaceutical, 32);
00414 scanStart=localtime(&img->scanStart);
00415 if(scanStart!=NULL) {
00416 main_header.scan_start_year=scanStart->tm_year+1900;
00417 main_header.scan_start_month=scanStart->tm_mon+1;
00418 main_header.scan_start_day=scanStart->tm_mday;
00419 main_header.scan_start_hour=scanStart->tm_hour;
00420 main_header.scan_start_minute=scanStart->tm_min;
00421 main_header.scan_start_second=scanStart->tm_sec;
00422 } else {
00423 main_header.scan_start_year=1900;
00424 main_header.scan_start_month=1;
00425 main_header.scan_start_day=1;
00426 main_header.scan_start_hour=0;
00427 main_header.scan_start_minute=0;
00428 main_header.scan_start_second=0;
00429 }
00430 main_header.isotope_halflife=img->isotopeHalflife;
00431 strcpy(main_header.isotope_code, imgIsotope(img));
00432 strcpy(main_header.study_name, img->studyNr);
00433 strcpy(main_header.patient_name, img->patientName);
00434 strcpy(main_header.patient_id, img->patientID);
00435 strncpy(main_header.user_process_code, img->userProcessCode, 10);
00436 strncpy(main_header.study_description, img->studyDescription, 32);
00437 if(IMG_TEST) ecat63PrintMainheader(&main_header, stdout);
00438
00439
00440 sdata=(short int*)malloc(img->dimx*img->dimy*sizeof(short int));
00441 if(sdata==NULL) {strcpy(ecat63errmsg, "out of memory"); return(4);}
00442
00443
00444 fp=ecat63Create(fname, &main_header);
00445 if(fp==NULL) {strcpy(ecat63errmsg, "cannot write ECAT file"); return(3);}
00446
00447
00448 switch(main_header.file_type) {
00449 case RAW_DATA:
00450 scan_header.data_type=main_header.data_type;
00451 scan_header.dimension_1=img->dimx;
00452 scan_header.dimension_2=img->dimy;
00453 scan_header.frame_duration_sec=1;
00454 scan_header.scale_factor=1.0;
00455 scan_header.frame_start_time=0;
00456 scan_header.frame_duration=1000;
00457 scan_header.loss_correction_fctr=1.0;
00458
00459 break;
00460 case IMAGE_DATA:
00461 image_header.data_type=main_header.data_type;
00462 image_header.num_dimensions=2;
00463 image_header.dimension_1=img->dimx;
00464 image_header.dimension_2=img->dimy;
00465 image_header.recon_scale=img->zoom;
00466 image_header.quant_scale=1.0;
00467 image_header.slice_width=img->sizez/10.;
00468 image_header.pixel_size=img->sizex/10.;
00469 image_header.frame_start_time=0;
00470 image_header.frame_duration=1000;
00471 image_header.plane_eff_corr_fctr=1.0;
00472 image_header.decay_corr_fctr=1.0;
00473 image_header.loss_corr_fctr=1.0;
00474 image_header.ecat_calibration_fctr=1.0;
00475 image_header.well_counter_cal_fctr=1.0;
00476 image_header.quant_units=main_header.calibration_units;
00477
00478 break;
00479 }
00480
00481
00482 n=0;
00483 for(plane=1; plane<=img->dimz; plane++) for(frame=1; frame<=img->dimt; frame++) {
00484
00485
00486 fmin=fmax=f=img->m[plane-1][0][0][frame-1];
00487 for(i=0; i<img->dimy; i++) for(j=0; j<img->dimx; j++) {
00488 f=img->m[plane-1][i][j][frame-1];
00489 if(f>fmax) fmax=f; else if(f<fmin) fmin=f;
00490 }
00491
00492 if(fabs(fmin)>fabs(fmax)) g=fabs(fmin); else g=fabs(fmax);
00493 if(g!=0) scale=32766./g; else scale=1.0;
00494
00495
00496 sptr=sdata;
00497 for(i=0; i<img->dimy; i++) for(j=0; j<img->dimx; j++) {
00498 *sptr=(short int)temp_roundf(scale*img->m[plane-1][i][j][frame-1]);
00499 sptr++;
00500 }
00501
00502 smin=(short int)temp_roundf(scale*fmin); smax=(short int)temp_roundf(scale*fmax);
00503 if(main_header.file_type==RAW_DATA) {
00504 scan_header.scan_min=smin; scan_header.scan_max=smax;
00505 scan_header.scale_factor=1.0/scale;
00506 scan_header.frame_start_time=(int)temp_roundf(1000.*img->start[frame-1]);
00507 scan_header.frame_duration=
00508 (int)temp_roundf(1000.*(img->end[frame-1]-img->start[frame-1]));
00509 scan_header.sample_distance=(img->sampleDistance)/10.0;
00510 scan_header.prompts=temp_roundf(img->prompts[frame-1]);
00511 scan_header.delayed=temp_roundf(img->randoms[frame-1]);
00512 } else if(main_header.file_type==IMAGE_DATA) {
00513 image_header.image_min=smin; image_header.image_max=smax;
00514 image_header.quant_scale=1.0/scale;
00515 image_header.frame_start_time=(int)temp_roundf(1000.*img->start[frame-1]);
00516 image_header.frame_duration=
00517 (int)temp_roundf(1000.*(img->end[frame-1]-img->start[frame-1]));
00518
00519 if(img->decayCorrected)
00520 image_header.decay_corr_fctr=img->decayCorrFactor[frame-1];
00521 }
00522
00523 m=mat_numcod(frame, img->planeNumber[plane-1], 1, 0, 0);
00524
00525 sptr=sdata;
00526 if(IMG_TEST) printf(" writing matnum=%d\n", m);
00527 if(main_header.file_type==RAW_DATA) {
00528 if(IMG_TEST) ecat63PrintScanheader(&scan_header, stdout);
00529 ret=ecat63WriteScan(fp, m, &scan_header, sptr);
00530 } else if(main_header.file_type==IMAGE_DATA) {
00531 if(IMG_TEST) ecat63PrintImageheader(&image_header, stdout);
00532 ret=ecat63WriteImage(fp, m, &image_header, sptr);
00533 }
00534 if(ret) {
00535 sprintf(ecat63errmsg, "cannot write data on pl%02d fr%02d (%d).",
00536 plane, frame, ret);
00537 fclose(fp); remove(fname); free(sdata);
00538 return(9);
00539 }
00540 n++;
00541 }
00542 if(IMG_TEST) printf(" %d matrices written in %s\n", n, fname);
00543
00544
00545 fclose(fp); free(sdata);
00546
00547 return(0);
00548 }
00549
00550
00551
00568 int ecat63ReadPlaneToImg(const char *fname, IMG *img) {
00569 int i, j, m, ret, blkNr=0, dim_x=0, dim_y=0, pxlNr=0, del_nr=0;
00570 int frameNr, planeNr, datatype=0, firstm=0, isFirst=0;
00571 int frame, plane, prev_frame, prev_plane, prev_frameNr=0;
00572 FILE *fp;
00573 ECAT63_mainheader main_header;
00574 MATRIXLIST mlist;
00575 MatDir tmpMatdir;
00576 Matval matval;
00577 ECAT63_imageheader image_header;
00578 ECAT63_scanheader scan_header;
00579 ECAT63_attnheader attn_header;
00580 ECAT63_normheader norm_header;
00581 float scale=1.0;
00582 short int *sptr;
00583 char *mdata=NULL, *mptr;
00584 int *iptr;
00585 struct tm scanStart;
00586
00587
00588 if(IMG_TEST) printf("ecat63ReadPlaneToImg(%s, *img)\n", fname);
00589
00590 if(fname==NULL) {strcpy(ecat63errmsg, "invalid ECAT filename"); return(2);}
00591 if(img==NULL || img->status==IMG_STATUS_UNINITIALIZED) {
00592 strcpy(ecat63errmsg, "image data not initialized"); return(2);}
00593
00594
00595 if((fp=fopen(fname, "rb")) == NULL) {
00596 strcpy(ecat63errmsg, "cannot open ECAT file"); return(3);}
00597
00598
00599 if((ret=ecat63ReadMainheader(fp, &main_header))) {
00600 sprintf(ecat63errmsg, "cannot read main header (%d)", ret);
00601 fclose(fp); return(4);
00602 }
00603 if(IMG_TEST) ecat63PrintMainheader(&main_header, stdout);
00604
00605
00606 ecat63InitMatlist(&mlist);
00607 ret=ecat63ReadMatlist(fp, &mlist);
00608 if(ret) {
00609 sprintf(ecat63errmsg, "cannot read matrix list (%d)", ret);
00610 fclose(fp); return(5);
00611 }
00612 if(mlist.matrixNr<=0) {
00613 strcpy(ecat63errmsg, "matrix list is empty"); fclose(fp); return(5);}
00614
00615 for(i=0; i<mlist.matrixNr-1; i++) for(j=i+1; j<mlist.matrixNr; j++) {
00616 if(mlist.matdir[i].matnum>mlist.matdir[j].matnum) {
00617 tmpMatdir=mlist.matdir[i];
00618 mlist.matdir[i]=mlist.matdir[j]; mlist.matdir[j]=tmpMatdir;
00619 }
00620 }
00621
00622 if(main_header.num_frames>0) {
00623 del_nr=ecat63DeleteLateFrames(&mlist, main_header.num_frames);
00624 if(IMG_TEST && del_nr>0)
00625 printf(" %d entries in matrix list will not be used.\n", del_nr);
00626 }
00627
00628
00629 if(img->status==IMG_STATUS_OCCUPIED) {
00630
00631 isFirst=0; prev_frameNr=img->dimt;
00632
00633 prev_plane=img->planeNumber[img->dimz-1];
00634
00635 } else {
00636
00637 isFirst=1;
00638
00639 prev_plane=-1;
00640 }
00641
00642 for(m=0, plane=-1; m<mlist.matrixNr; m++) if(mlist.matdir[m].matstat==1) {
00643 mat_numdoc(mlist.matdir[m].matnum, &matval);
00644 if(matval.plane>prev_plane) {plane=matval.plane; break;}
00645 }
00646
00647 if(plane<0) {
00648 fclose(fp); ecat63EmptyMatlist(&mlist);
00649 if(isFirst) {
00650 sprintf(ecat63errmsg, "ECAT file contains no matrices");
00651 return(7);
00652 } else {
00653 sprintf(ecat63errmsg, "ECAT file contains no more planes");
00654 if(IMG_TEST) printf("%s\n", ecat63errmsg);
00655 return(1);
00656 }
00657 }
00658 if(IMG_TEST) printf("Next plane to read: %d\n", plane);
00659
00660 imgEmpty(img);
00661
00662
00663 prev_frame=frame=-1; frameNr=0; ret=0; planeNr=1;
00664 for(m=0; m<mlist.matrixNr; m++) if(mlist.matdir[m].matstat==1) {
00665
00666 mat_numdoc(mlist.matdir[m].matnum, &matval);
00667 if(matval.plane<plane) continue; else if(matval.plane>plane) break;
00668 if(main_header.num_frames>=main_header.num_gates) frame=matval.frame;
00669 else frame=matval.gate;
00670 frameNr++;
00671
00672 if(prev_frame>0 && frame!=prev_frame+1) {ret=1; break;}
00673 prev_frame=frame;
00674
00675 if(frameNr==1) {
00676 firstm=m;
00677 blkNr=mlist.matdir[m].endblk-mlist.matdir[m].strtblk;
00678 } else if(blkNr!=mlist.matdir[m].endblk-mlist.matdir[m].strtblk) {
00679 ret=2; break;
00680 }
00681 prev_frame=frame;
00682 }
00683 if(ret==1) {
00684 strcpy(ecat63errmsg, "missing matrix");
00685 ecat63EmptyMatlist(&mlist); fclose(fp); return(6);
00686 }
00687 if(ret==2) {
00688 strcpy(ecat63errmsg, "matrix sizes are different");
00689 ecat63EmptyMatlist(&mlist); fclose(fp); return(6);
00690 }
00691
00692 if(!isFirst && frameNr!=prev_frameNr) {
00693 strcpy(ecat63errmsg, "different frame nr in different planes");
00694 ecat63EmptyMatlist(&mlist); fclose(fp); return(6);
00695 }
00696
00697
00698 m=firstm; if(main_header.file_type==IMAGE_DATA) {
00699 ret=ecat63ReadImageheader(fp, mlist.matdir[m].strtblk, &image_header);
00700 dim_x=image_header.dimension_1; dim_y=image_header.dimension_2;
00701 } else if(main_header.file_type==RAW_DATA) {
00702 ret=ecat63ReadScanheader(fp, mlist.matdir[m].strtblk, &scan_header);
00703 dim_x=scan_header.dimension_1; dim_y=scan_header.dimension_2;
00704 } else if(main_header.file_type==ATTN_DATA) {
00705 ret=ecat63ReadAttnheader(fp, mlist.matdir[m].strtblk, &attn_header);
00706 dim_x=attn_header.dimension_1; dim_y=attn_header.dimension_2;
00707 } else if(main_header.file_type==NORM_DATA) {
00708 ret=ecat63ReadNormheader(fp, mlist.matdir[m].strtblk, &norm_header);
00709 dim_x=norm_header.dimension_1; dim_y=norm_header.dimension_2;
00710 }
00711 pxlNr=dim_x*dim_y;
00712 if(ret) {
00713 sprintf(ecat63errmsg, "cannot read matrix %u subheader", mlist.matdir[m].matnum);
00714 ecat63EmptyMatlist(&mlist); fclose(fp); return(7);
00715 }
00716
00717
00718 if(IMG_TEST) printf("allocating memory for %d blocks\n", blkNr);
00719 mdata=(char*)malloc(blkNr*MatBLKSIZE); if(mdata==NULL) {
00720 strcpy(ecat63errmsg, "out of memory");
00721 fclose(fp); ecat63EmptyMatlist(&mlist); return(8);
00722 }
00723
00724 ret=imgAllocate(img, planeNr, dim_y, dim_x, frameNr);
00725 if(ret) {
00726 sprintf(ecat63errmsg, "out of memory (%d)", ret);
00727 fclose(fp); ecat63EmptyMatlist(&mlist); free(mdata); return(8);
00728 }
00729
00730
00731 img->scanner=main_header.system_type;
00732 img->unit=main_header.calibration_units;
00733 strncpy(img->radiopharmaceutical, main_header.radiopharmaceutical, 32);
00734 img->isotopeHalflife=main_header.isotope_halflife;
00735 memset(&scanStart, 0, sizeof(struct tm));
00736 scanStart.tm_year=main_header.scan_start_year-1900;
00737 scanStart.tm_mon=main_header.scan_start_month-1;
00738 scanStart.tm_mday=main_header.scan_start_day; scanStart.tm_yday=0;
00739 scanStart.tm_hour=main_header.scan_start_hour;
00740 scanStart.tm_min=main_header.scan_start_minute;
00741 scanStart.tm_sec=main_header.scan_start_second;
00742 scanStart.tm_isdst=-1;
00743 img->scanStart=mktime(&scanStart);
00744 if(img->scanStart==-1) img->scanStart=0;
00745 img->axialFOV=10.*main_header.axial_fov;
00746 img->transaxialFOV=10.*main_header.transaxial_fov;
00747 strncpy(img->studyNr, main_header.study_name, MAX_STUDYNR_LEN);
00748 strncpy(img->patientName, main_header.patient_name, 31);
00749 strncpy(img->patientID, main_header.patient_id, 15);
00750 img->sizez=10.*main_header.plane_separation;
00751 if(main_header.file_type==IMAGE_DATA) {
00752 img->type=IMG_TYPE_IMAGE;
00753 if(img->unit<1) img->unit=image_header.quant_units;
00754 img->zoom=image_header.recon_scale;
00755 if(image_header.decay_corr_fctr>1.0) img->decayCorrected=1;
00756 img->sizex=img->sizey=10.*image_header.pixel_size;
00757 if(img->sizez<10.*image_header.slice_width)
00758 img->sizez=10.*image_header.slice_width;
00759 } else if(main_header.file_type==RAW_DATA) {
00760 img->type=IMG_TYPE_RAW;
00761 img->decayCorrected=0;
00762 } else if(main_header.file_type==ATTN_DATA) {
00763 img->type=IMG_TYPE_RAW;
00764 img->decayCorrected=0;
00765 } else if(main_header.file_type==NORM_DATA) {
00766 img->type=IMG_TYPE_RAW;
00767 img->decayCorrected=0;
00768 }
00769 img->planeNumber[0]=plane;
00770 strncpy(img->studyDescription, main_header.study_description, 32);
00771 strncpy(img->userProcessCode, main_header.user_process_code, 10); img->userProcessCode[10]=(char)0;
00772
00773 if(!img->studyNr[0] && studynr_validity_check(img->userProcessCode))
00774 strcpy(img->studyNr, img->userProcessCode);
00775
00776 img->_fileFormat=IMG_E63;
00777
00778
00779 for(m=firstm, frame=1; m<mlist.matrixNr; m++) if(mlist.matdir[m].matstat==1) {
00780
00781 mat_numdoc(mlist.matdir[m].matnum, &matval);
00782 if(matval.plane>plane) break;
00783
00784 if(main_header.file_type==IMAGE_DATA) {
00785 ret=ecat63ReadImageheader(fp, mlist.matdir[m].strtblk, &image_header);
00786 if(dim_x!=image_header.dimension_1 || dim_y!=image_header.dimension_2) ret=1;
00787 scale=image_header.quant_scale;
00788 if(image_header.ecat_calibration_fctr>0.0
00789 && fabs(main_header.calibration_factor/image_header.ecat_calibration_fctr-1.0)>0.0001)
00790 scale*=image_header.ecat_calibration_fctr;
00791 if(img->unit==0 && image_header.quant_units>0) img->unit=image_header.quant_units;
00792 datatype=image_header.data_type;
00793 img->start[frame-1]=image_header.frame_start_time/1000.;
00794 img->end[frame-1]=img->start[frame-1]+image_header.frame_duration/1000.;
00795 img->mid[frame-1]=0.5*(img->start[frame-1]+img->end[frame-1]);
00796 if(image_header.decay_corr_fctr>1.0)
00797 img->decayCorrFactor[frame-1]=image_header.decay_corr_fctr;
00798 } else if(main_header.file_type==RAW_DATA) {
00799 ret=ecat63ReadScanheader(fp, mlist.matdir[m].strtblk, &scan_header);
00800 if(dim_x!=scan_header.dimension_1 || dim_y!=scan_header.dimension_2) ret=1;
00801 scale=scan_header.scale_factor;
00802 if(scan_header.loss_correction_fctr>0.0) scale*=scan_header.loss_correction_fctr;
00803 datatype=scan_header.data_type;
00804 img->start[frame-1]=scan_header.frame_start_time/1000.;
00805 img->end[frame-1]=img->start[frame-1]+scan_header.frame_duration/1000.;
00806 img->mid[frame-1]=0.5*(img->start[frame-1]+img->end[frame-1]);
00807 img->sampleDistance=10.0*scan_header.sample_distance;
00808 img->prompts[frame-1]=(float)scan_header.prompts;
00809 img->randoms[frame-1]=(float)scan_header.delayed;
00810 } else if(main_header.file_type==ATTN_DATA) {
00811 ret=ecat63ReadAttnheader(fp, mlist.matdir[m].strtblk, &attn_header);
00812 if(dim_x!=attn_header.dimension_1 || dim_y!=attn_header.dimension_2) ret=1;
00813 scale=attn_header.scale_factor;
00814 datatype=attn_header.data_type;
00815 } else if(main_header.file_type==NORM_DATA) {
00816 ret=ecat63ReadNormheader(fp, mlist.matdir[m].strtblk, &norm_header);
00817 if(dim_x!=norm_header.dimension_1 || dim_y!=norm_header.dimension_2) ret=1;
00818 scale=norm_header.scale_factor;
00819 datatype=norm_header.data_type;
00820 } else datatype=-1;
00821 if(ret) {
00822 sprintf(ecat63errmsg, "cannot read matrix %u subheader", mlist.matdir[m].matnum);
00823 ecat63EmptyMatlist(&mlist); free(mdata); fclose(fp); return(7);
00824 }
00825 if(main_header.calibration_factor>0.0) scale*=main_header.calibration_factor;
00826 if(IMG_TEST>2) printf("scale=%e datatype=%d\n", scale, datatype);
00827
00828 ret=ecat63ReadMatdata(fp, mlist.matdir[m].strtblk+1,
00829 mlist.matdir[m].endblk-mlist.matdir[m].strtblk,
00830 mdata, datatype);
00831 if(ret) {
00832 strcpy(ecat63errmsg, "cannot read matrix data");
00833 ecat63EmptyMatlist(&mlist); free(mdata); fclose(fp); return(9);
00834 }
00835 if(datatype==BYTE_TYPE) {
00836 for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr++) {
00837 img->m[0][i][j][frame-1]=scale*(float)(*mptr);
00838 }
00839 } else if(datatype==VAX_I2 || datatype==SUN_I2) {
00840 for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr+=2) {
00841 sptr=(short int*)mptr;
00842 img->m[0][i][j][frame-1]=scale*(float)(*sptr);
00843 }
00844 } else if(datatype==VAX_I4 || datatype==SUN_I4) {
00845 for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr+=4) {
00846 iptr=(int*)mptr;
00847 img->m[0][i][j][frame-1]=scale*(float)(*iptr);
00848 }
00849 } else if(datatype==VAX_R4 || datatype==IEEE_R4) {
00850 for(i=0, mptr=mdata; i<dim_y; i++) for(j=0; j<dim_x; j++, mptr+=4) {
00851 memcpy(&img->m[0][i][j][frame-1], mptr, 4);
00852 img->m[0][i][j][frame-1]*=scale;
00853 }
00854 }
00855 frame++;
00856 }
00857
00858 imgUnitFromEcat(img, img->unit);
00859
00860
00861 free(mdata);
00862 ecat63EmptyMatlist(&mlist);
00863 fclose(fp);
00864
00865
00866 if(img->_dataType<1) img->_dataType=datatype;
00867
00868 if(img->_dataType==VAX_I4 || img->_dataType==VAX_R4) img->_dataType=VAX_I2;
00869
00870 return(0);
00871 }
00872
00873
00874
00886 int ecat63AddImg(const char *fname, IMG *img) {
00887 int n, i, j, m, ret=0, add=0;
00888 int frameNr, planeNr;
00889 int frame, plane, prev_frame, prev_plane;
00890 float f, fmax, fmin, g, scale;
00891 short int *sdata, *sptr, smin, smax;
00892 FILE *fp;
00893 ECAT63_mainheader main_header;
00894 ECAT63_imageheader image_header;
00895 ECAT63_scanheader scan_header;
00896 MATRIXLIST mlist;
00897 MatDir tmpMatdir;
00898 Matval matval;
00899 struct tm *scanStart;
00900
00901
00902 if(IMG_TEST) printf("ecat63AddImg(%s, *img)\n", fname);
00903 if(IMG_TEST>1) imgInfo(img);
00904
00905 if(fname==NULL) {strcpy(ecat63errmsg, "invalid ECAT filename"); return(1);}
00906 if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) {
00907 strcpy(ecat63errmsg, "image data is empty"); return(2);}
00908 if(img->_dataType<1) img->_dataType=VAX_I2;
00909
00910
00911 memset(&main_header, 0, sizeof(ECAT63_mainheader));
00912 memset(&image_header,0, sizeof(ECAT63_imageheader));
00913 memset(&scan_header, 0, sizeof(ECAT63_scanheader));
00914
00915
00916 main_header.sw_version=2;
00917 main_header.num_planes=img->dimz;
00918 main_header.num_frames=img->dimt;
00919 main_header.num_gates=1;
00920 main_header.num_bed_pos=1;
00921 if(img->type==IMG_TYPE_IMAGE) main_header.file_type=IMAGE_DATA;
00922 else main_header.file_type=RAW_DATA;
00923 main_header.data_type=img->_dataType;
00924 if(img->scanner>0) main_header.system_type=img->scanner;
00925 else main_header.system_type=ECAT63_SYSTEM_TYPE_DEFAULT;
00926 main_header.calibration_factor=1.0;
00927 main_header.axial_fov=img->axialFOV/10.0;
00928 main_header.transaxial_fov=img->transaxialFOV/10.0;
00929 main_header.plane_separation=img->sizez/10.0;
00930 main_header.calibration_units=imgUnitToEcat6(img);
00931 strncpy(main_header.radiopharmaceutical, img->radiopharmaceutical, 32);
00932 scanStart=localtime(&img->scanStart);
00933 if(scanStart!=NULL) {
00934 main_header.scan_start_year=scanStart->tm_year+1900;
00935 main_header.scan_start_month=scanStart->tm_mon+1;
00936 main_header.scan_start_day=scanStart->tm_mday;
00937 main_header.scan_start_hour=scanStart->tm_hour;
00938 main_header.scan_start_minute=scanStart->tm_min;
00939 main_header.scan_start_second=scanStart->tm_sec;
00940 } else {
00941 main_header.scan_start_year=1900;
00942 main_header.scan_start_month=1;
00943 main_header.scan_start_day=1;
00944 main_header.scan_start_hour=0;
00945 main_header.scan_start_minute=0;
00946 main_header.scan_start_second=0;
00947 }
00948 main_header.isotope_halflife=img->isotopeHalflife;
00949 strcpy(main_header.isotope_code, imgIsotope(img));
00950 strcpy(main_header.study_name, img->studyNr);
00951 strcpy(main_header.patient_name, img->patientName);
00952 strcpy(main_header.patient_id, img->patientID);
00953 strncpy(main_header.study_description, img->studyDescription, 32);
00954 strncpy(main_header.user_process_code, img->userProcessCode, 10);
00955
00956
00957 if(access(fname, 0) != -1) {
00958 if(IMG_TEST) printf("Opening existing ECAT file.\n");
00959 add=1;
00960
00961 if((fp=fopen(fname, "rb+")) == NULL) {
00962 strcpy(ecat63errmsg, "cannot create ECAT file"); return(3);}
00963
00964 if((ret=ecat63ReadMainheader(fp, &main_header))) {
00965 sprintf(ecat63errmsg, "cannot read main header (%d)", ret);
00966 fclose(fp); return(3);
00967 }
00968 fflush(fp);
00969
00970 if((main_header.file_type==IMAGE_DATA && img->type==IMG_TYPE_IMAGE) ||
00971 (main_header.file_type==RAW_DATA && img->type==IMG_TYPE_RAW)) {
00972
00973 } else {
00974 sprintf(ecat63errmsg, "cannot add different filetype");
00975 fclose(fp); return(3);
00976 }
00977 } else {
00978 if(IMG_TEST) printf("ECAT file does not exist.\n");
00979 add=0;
00980
00981 fp=ecat63Create(fname, &main_header);
00982 if(fp==NULL) {strcpy(ecat63errmsg, "cannot create ECAT file"); return(3);}
00983 }
00984 if(IMG_TEST) ecat63PrintMainheader(&main_header, stdout);
00985
00986
00987 sdata=(short int*)malloc(img->dimx*img->dimy*sizeof(short int));
00988 if(sdata==NULL) {strcpy(ecat63errmsg, "out of memory"); return(4);}
00989
00990
00991 switch(main_header.file_type) {
00992 case RAW_DATA:
00993 scan_header.data_type=main_header.data_type;
00994 scan_header.dimension_1=img->dimx;
00995 scan_header.dimension_2=img->dimy;
00996 scan_header.frame_duration_sec=1;
00997 scan_header.scale_factor=1.0;
00998 scan_header.frame_start_time=0;
00999 scan_header.frame_duration=1000;
01000 scan_header.loss_correction_fctr=1.0;
01001 scan_header.sample_distance=(img->sampleDistance)/10.0;
01002
01003 break;
01004 case IMAGE_DATA:
01005 image_header.data_type=main_header.data_type;
01006 image_header.num_dimensions=2;
01007 image_header.dimension_1=img->dimx;
01008 image_header.dimension_2=img->dimy;
01009 image_header.recon_scale=img->zoom;
01010 image_header.quant_scale=1.0;
01011 image_header.slice_width=img->sizez/10.;
01012 image_header.pixel_size=img->sizex/10.;
01013 image_header.frame_start_time=0;
01014 image_header.frame_duration=1000;
01015 image_header.plane_eff_corr_fctr=1.0;
01016 image_header.decay_corr_fctr=0.0;
01017 image_header.loss_corr_fctr=1.0;
01018 image_header.ecat_calibration_fctr=1.0;
01019 image_header.well_counter_cal_fctr=1.0;
01020 image_header.quant_units=main_header.calibration_units;
01021
01022 break;
01023 }
01024
01025
01026 n=0;
01027 for(plane=1; plane<=img->dimz; plane++) for(frame=1; frame<=img->dimt; frame++) {
01028
01029
01030 fmin=fmax=f=img->m[plane-1][0][0][frame-1];
01031 for(i=0; i<img->dimy; i++) for(j=0; j<img->dimx; j++) {
01032 f=img->m[plane-1][i][j][frame-1];
01033 if(f>fmax) fmax=f; else if(f<fmin) fmin=f;
01034 }
01035
01036 if(fabs(fmin)>fabs(fmax)) g=fabs(fmin); else g=fabs(fmax);
01037 if(g!=0) scale=32766./g; else scale=1.0;
01038
01039
01040 sptr=sdata;
01041 for(i=0; i<img->dimy; i++) for(j=0; j<img->dimx; j++) {
01042 *sptr=(short int)temp_roundf(scale*img->m[plane-1][i][j][frame-1]);
01043 sptr++;
01044 }
01045
01046 smin=(short int)temp_roundf(scale*fmin); smax=(short int)temp_roundf(scale*fmax);
01047 if(main_header.file_type==RAW_DATA) {
01048 scan_header.scan_min=smin; scan_header.scan_max=smax;
01049 scan_header.scale_factor=1.0/scale;
01050 scan_header.frame_start_time=(int)temp_roundf(1000.*img->start[frame-1]);
01051 scan_header.frame_duration=
01052 (int)temp_roundf(1000.*(img->end[frame-1]-img->start[frame-1]));
01053 scan_header.prompts=temp_roundf(img->prompts[frame-1]);
01054 scan_header.delayed=temp_roundf(img->randoms[frame-1]);
01055 } else if(main_header.file_type==IMAGE_DATA) {
01056 image_header.image_min=smin; image_header.image_max=smax;
01057 image_header.quant_scale=1.0/scale;
01058 image_header.frame_start_time=(int)temp_roundf(1000.*img->start[frame-1]);
01059 image_header.frame_duration=
01060 (int)temp_roundf(1000.*(img->end[frame-1]-img->start[frame-1]));
01061
01062 if(img->decayCorrected)
01063 image_header.decay_corr_fctr=img->decayCorrFactor[frame-1];
01064 }
01065
01066 m=mat_numcod(frame, img->planeNumber[plane-1], 1, 0, 0);
01067 sptr=sdata;
01068 if(IMG_TEST) printf(" writing matnum=%d\n", m);
01069 if(main_header.file_type==RAW_DATA) {
01070
01071 ret=ecat63WriteScan(fp, m, &scan_header, sptr);
01072 } else if(main_header.file_type==IMAGE_DATA) {
01073
01074 ret=ecat63WriteImage(fp, m, &image_header, sptr);
01075 }
01076 if(ret) {
01077 sprintf(ecat63errmsg, "cannot write data on pl%02d fr%02d (%d).",
01078 plane, frame, ret);
01079 fclose(fp); remove(fname); free(sdata);
01080 return(9);
01081 }
01082 n++;
01083 }
01084 if(IMG_TEST) printf(" %d matrices written in %s\n", n, fname);
01085
01086
01087 free(sdata);
01088
01089
01090 if(add==1) {
01091 fflush(fp);
01092
01093 ecat63InitMatlist(&mlist);
01094 ret=ecat63ReadMatlist(fp, &mlist);
01095 if(ret) {
01096 sprintf(ecat63errmsg, "cannot read matrix list (%d)", ret);
01097 fclose(fp); return(21);
01098 }
01099 if(mlist.matrixNr<=0) {
01100 strcpy(ecat63errmsg, "matrix list is empty"); fclose(fp); return(21);}
01101
01102 for(i=0; i<mlist.matrixNr-1; i++) for(j=i+1; j<mlist.matrixNr; j++) {
01103 if(mlist.matdir[i].matnum>mlist.matdir[j].matnum) {
01104 tmpMatdir=mlist.matdir[i];
01105 mlist.matdir[i]=mlist.matdir[j]; mlist.matdir[j]=tmpMatdir;
01106 }
01107 }
01108
01109 prev_plane=plane=-1; prev_frame=frame=-1; frameNr=planeNr=0;
01110 for(m=0; m<mlist.matrixNr; m++) {
01111 mat_numdoc(mlist.matdir[m].matnum, &matval);
01112 plane=matval.plane; frame=matval.frame;
01113 if(plane!=prev_plane) {frameNr=1; planeNr++;} else {frameNr++;}
01114 prev_plane=plane; prev_frame=frame;
01115 }
01116 ecat63EmptyMatlist(&mlist);
01117 main_header.num_planes=planeNr; main_header.num_frames=frameNr;
01118
01119 ret=ecat63WriteMainheader(fp, &main_header);
01120 if(ret) {
01121 sprintf(ecat63errmsg, "cannot write mainheader (%d)", ret);
01122 fclose(fp); return(22);
01123 }
01124 }
01125
01126
01127 fclose(fp);
01128
01129 return(0);
01130 }
01131
01132
01133
01140 int imgEcat63Supported(ECAT63_mainheader *h) {
01141 if(h==NULL) return(0);
01142 if(h->file_type==IMAGE_DATA) return(1);
01143 if(h->file_type==RAW_DATA) return(1);
01144 if(h->file_type==ATTN_DATA) return(1);
01145 if(h->file_type==NORM_DATA) return(1);
01146 return(0);
01147 }
01148
01149
01150
01157 void imgGetEcat63MHeader(IMG *img, ECAT63_mainheader *h) {
01158 struct tm scanStart;
01159
01160 img->_dataType=h->data_type;
01161
01162 if(img->_dataType==VAX_I4 || img->_dataType==VAX_R4) img->_dataType=VAX_I2;
01163 img->scanner=h->system_type;
01164 strncpy(img->radiopharmaceutical, h->radiopharmaceutical, 32);
01165 img->isotopeHalflife=h->isotope_halflife;
01166 {
01167 memset(&scanStart, 0, sizeof(struct tm));
01168 scanStart.tm_year=h->scan_start_year-1900;
01169 scanStart.tm_mon=h->scan_start_month-1;
01170 scanStart.tm_mday=h->scan_start_day; scanStart.tm_yday=0;
01171 scanStart.tm_hour=h->scan_start_hour;
01172 scanStart.tm_min=h->scan_start_minute;
01173 scanStart.tm_sec=h->scan_start_second;
01174 scanStart.tm_isdst=-1;
01175 img->scanStart=mktime(&scanStart); if(img->scanStart==-1) img->scanStart=0;
01176 }
01177 img->axialFOV=10.0*h->axial_fov;
01178 img->transaxialFOV=10.0*h->transaxial_fov;
01179 strncpy(img->studyNr, h->study_name, MAX_STUDYNR_LEN);
01180 strncpy(img->patientName, h->patient_name, 31);
01181 strncpy(img->patientID, h->patient_id, 15);
01182 img->sizez=10.0*h->plane_separation;
01183 switch(h->file_type) {
01184 case IMAGE_DATA: img->type=IMG_TYPE_IMAGE; break;
01185 case RAW_DATA:
01186 case ATTN_DATA:
01187 case NORM_DATA:
01188 default: img->type=IMG_TYPE_RAW;
01189 }
01190 strncpy(img->studyDescription, h->study_description, 32);
01191 strncpy(img->userProcessCode, h->user_process_code, 10);
01192 img->userProcessCode[10]=(char)0;
01193
01194 if(!img->studyNr[0] && studynr_validity_check(img->userProcessCode))
01195 strcpy(img->studyNr, img->userProcessCode);
01196
01197 imgUnitFromEcat(img, h->calibration_units);
01198 }
01199
01200
01201
01208 void imgSetEcat63MHeader(IMG *img, ECAT63_mainheader *h) {
01209 struct tm *scanStart;
01210
01211 h->sw_version=2;
01212 h->num_planes=img->dimz;
01213 h->num_frames=img->dimt;
01214 h->num_gates=1;
01215 h->num_bed_pos=1;
01216 if(img->type==IMG_TYPE_IMAGE) h->file_type=IMAGE_DATA;
01217 else h->file_type=RAW_DATA;
01218 h->data_type=VAX_I2;
01219 if(img->scanner>0) h->system_type=img->scanner;
01220 else h->system_type=ECAT63_SYSTEM_TYPE_DEFAULT;
01221 h->calibration_factor=1.0;
01222 h->axial_fov=img->axialFOV/10.0;
01223 h->transaxial_fov=img->transaxialFOV/10.0;
01224 h->plane_separation=img->sizez/10.0;
01225 h->calibration_units=imgUnitToEcat6(img);
01226 strncpy(h->radiopharmaceutical, img->radiopharmaceutical, 32);
01227 scanStart=localtime(&img->scanStart);
01228 if(scanStart!=NULL) {
01229 h->scan_start_year=scanStart->tm_year+1900;
01230 h->scan_start_month=scanStart->tm_mon+1;
01231 h->scan_start_day=scanStart->tm_mday;
01232 h->scan_start_hour=scanStart->tm_hour;
01233 h->scan_start_minute=scanStart->tm_min;
01234 h->scan_start_second=scanStart->tm_sec;
01235 } else {
01236 h->scan_start_year=1900;
01237 h->scan_start_month=1;
01238 h->scan_start_day=1;
01239 h->scan_start_hour=0;
01240 h->scan_start_minute=0;
01241 h->scan_start_second=0;
01242 }
01243 h->isotope_halflife=img->isotopeHalflife;
01244 strcpy(h->isotope_code, imgIsotope(img));
01245 strcpy(h->study_name, img->studyNr);
01246 strcpy(h->patient_name, img->patientName);
01247 strcpy(h->patient_id, img->patientID);
01248 strncpy(h->user_process_code, img->userProcessCode, 10);
01249 strncpy(h->study_description, img->studyDescription, 32);
01250 }
01251
01252
01253
01260 int imgGetEcat63Fileformat(ECAT63_mainheader *h) {
01261 int fileFormat=IMG_UNKNOWN;
01262 switch(h->file_type) {
01263 case IMAGE_DATA:
01264 case RAW_DATA:
01265 case ATTN_DATA:
01266 case NORM_DATA:
01267 fileFormat=IMG_E63; break;
01268 default:
01269 fileFormat=IMG_UNKNOWN; break;
01270 }
01271 return fileFormat;
01272 }
01273
01274
01275
01289 int imgReadEcat63Header(const char *fname, IMG *img) {
01290 int m, blkNr=0, ret=0;
01291 int frameNr, planeNr, del_nr=0;
01292 FILE *fp;
01293 ECAT63_mainheader main_header;
01294 MATRIXLIST mlist;
01295 ECAT63_imageheader image_header;
01296 ECAT63_scanheader scan_header;
01297 ECAT63_attnheader attn_header;
01298 ECAT63_normheader norm_header;
01299
01300
01301 if(IMG_TEST) printf("\nimgReadEcat63Header(%s, *img)\n", fname);
01302
01303
01304 if(img==NULL) return STATUS_FAULT;
01305 if(img->status!=IMG_STATUS_INITIALIZED) return STATUS_FAULT;
01306 imgSetStatus(img, STATUS_FAULT);
01307 if(fname==NULL) return STATUS_FAULT;
01308
01309
01310 if((fp=fopen(fname, "rb")) == NULL) return STATUS_NOFILE;
01311
01312
01313 if((ret=ecat63ReadMainheader(fp, &main_header))) {
01314 fclose(fp); return STATUS_NOMAINHEADER;}
01315
01316 if(imgEcat63Supported(&main_header)==0) {fclose(fp); return STATUS_UNSUPPORTED;}
01317
01318 imgGetEcat63MHeader(img, &main_header);
01319 if(IMG_TEST>7) printf("img.type := %d\n", img->type);
01320 img->_fileFormat=imgGetEcat63Fileformat(&main_header);
01321 if(IMG_TEST>7) printf("img._fileFormat := %d\n", img->_fileFormat);
01322 if(img->_fileFormat==IMG_UNKNOWN) {fclose(fp); return STATUS_UNSUPPORTED;}
01323
01324
01325 ecat63InitMatlist(&mlist);
01326 ret=ecat63ReadMatlist(fp, &mlist);
01327 if(ret) {fclose(fp); return STATUS_NOMATLIST;}
01328 if(mlist.matrixNr<=0) {
01329 ecat63EmptyMatlist(&mlist); fclose(fp); return STATUS_INVALIDMATLIST;}
01330
01331 ecat63GatherMatlist(&mlist, 1, 1, 1, 1);
01332 ecat63SortMatlistByPlane(&mlist);
01333
01334 if(main_header.num_frames>0)
01335 del_nr=ecat63DeleteLateFrames(&mlist, main_header.num_frames);
01336
01337 ret=ecat63GetPlaneAndFrameNr(&mlist, &main_header, &planeNr, &frameNr);
01338 if(ret) {ecat63EmptyMatlist(&mlist); fclose(fp); return ret;}
01339 img->dimz=planeNr;
01340 img->dimt=frameNr;
01341
01342 ret=ecat63GetMatrixBlockSize(&mlist, &blkNr);
01343 if(ret) {ecat63EmptyMatlist(&mlist); fclose(fp); return ret;}
01344
01345
01346 if(IMG_TEST>5) printf("main_header.file_type := %d\n", main_header.file_type);
01347 m=0;
01348 switch(main_header.file_type) {
01349 case IMAGE_DATA:
01350 ret=ecat63ReadImageheader(fp, mlist.matdir[m].strtblk, &image_header);
01351 break;
01352 case RAW_DATA:
01353 ret=ecat63ReadScanheader(fp, mlist.matdir[m].strtblk, &scan_header);
01354 break;
01355 case ATTN_DATA:
01356 ret=ecat63ReadAttnheader(fp, mlist.matdir[m].strtblk, &attn_header);
01357 break;
01358 case NORM_DATA:
01359 ret=ecat63ReadNormheader(fp, mlist.matdir[m].strtblk, &norm_header);
01360 break;
01361 default: ret=-1;
01362 }
01363
01364 ecat63EmptyMatlist(&mlist); fclose(fp);
01365
01366 if(ret) return STATUS_NOSUBHEADER;
01367
01368
01369
01370
01371
01372
01373 switch(main_header.file_type) {
01374 case IMAGE_DATA:
01375 img->dimx=image_header.dimension_1; img->dimy=image_header.dimension_2;
01376 if(img->unit<1) imgUnitFromEcat(img, image_header.quant_units);
01377 img->_dataType=image_header.data_type;
01378 img->zoom=image_header.recon_scale;
01379 if(image_header.decay_corr_fctr>1.0) img->decayCorrected=1;
01380 img->sizex=img->sizey=10.*image_header.pixel_size;
01381 if(img->sizez<10.*image_header.slice_width)
01382 img->sizez=10.*image_header.slice_width;
01383 break;
01384 case RAW_DATA:
01385 img->dimx=scan_header.dimension_1; img->dimy=scan_header.dimension_2;
01386 img->type=IMG_TYPE_RAW;
01387 img->_dataType=scan_header.data_type;
01388 img->decayCorrected=0;
01389 img->sampleDistance=10.0*scan_header.sample_distance;
01390 break;
01391 case ATTN_DATA:
01392 img->dimx=attn_header.dimension_1; img->dimy=attn_header.dimension_2;
01393 img->type=IMG_TYPE_RAW;
01394 img->decayCorrected=0;
01395 img->_dataType=attn_header.data_type;
01396 break;
01397 case NORM_DATA:
01398 img->dimx=norm_header.dimension_1; img->dimy=norm_header.dimension_2;
01399 img->type=IMG_TYPE_RAW;
01400 img->decayCorrected=0;
01401 img->_dataType=norm_header.data_type;
01402 break;
01403 }
01404
01405
01406 if(img->_dataType==VAX_I4 || img->_dataType==VAX_R4) img->_dataType=VAX_I2;
01407
01408 imgSetStatus(img, STATUS_OK);
01409 return STATUS_OK;
01410 }
01411
01412
01413
01422 int imgReadEcat63FirstFrame(const char *fname, IMG *img) {
01423 int ret=0;
01424
01425 if(IMG_TEST) printf("\nimgReadEcat63FirstFrame(%s, *img)\n", fname);
01426
01427 if(img==NULL) return STATUS_FAULT;
01428 if(img->status!=IMG_STATUS_INITIALIZED) return STATUS_FAULT;
01429 imgSetStatus(img, STATUS_FAULT);
01430 if(fname==NULL) return STATUS_FAULT;
01431
01432
01433 ret=imgReadEcat63Header(fname, img); if(ret) return(ret);
01434 if(IMG_TEST>3) imgInfo(img);
01435
01436
01437 img->dimt=1;
01438 ret=imgAllocate(img, img->dimz, img->dimy, img->dimx, img->dimt);
01439 if(ret) return STATUS_NOMEMORY;
01440
01441
01442 ret=imgReadEcat63Frame(fname, 1, img, 0); if(ret) return(ret);
01443
01444
01445 imgSetStatus(img, STATUS_OK);
01446 return STATUS_OK;
01447 }
01448
01449
01450
01464 int imgReadEcat63Frame(const char *fname, int frame_to_read, IMG *img, int frame_index) {
01465 int m, ret=0, blkNr=0, frame, plane, seqplane=-1, xi, yi;
01466 int local_data_type=0;
01467 FILE *fp;
01468 ECAT63_mainheader main_header;
01469 MATRIXLIST mlist;
01470 Matval matval;
01471 ECAT63_imageheader image_header;
01472 ECAT63_scanheader scan_header;
01473 ECAT63_attnheader attn_header;
01474 ECAT63_normheader norm_header;
01475 float scale=1.0;
01476 short int *sptr;
01477 char *mdata=NULL, *mptr;
01478 int *iptr;
01479
01480
01481 if(IMG_TEST) printf("\nimgReadEcat63Frame(%s, %d, *img, %d)\n",
01482 fname, frame_to_read, frame_index);
01483
01484
01485 if(img==NULL) return STATUS_FAULT;
01486 if(img->status!=IMG_STATUS_OCCUPIED) return STATUS_FAULT;
01487 imgSetStatus(img, STATUS_FAULT);
01488 if(fname==NULL) return STATUS_FAULT;
01489 if(frame_index<0 || frame_index>img->dimt-1) return STATUS_FAULT;
01490 if(frame_to_read<1) return STATUS_FAULT;
01491
01492
01493 if((fp=fopen(fname, "rb")) == NULL) return STATUS_NOFILE;
01494
01495
01496 if((ret=ecat63ReadMainheader(fp, &main_header))) {
01497 fclose(fp); return STATUS_NOMAINHEADER;}
01498
01499
01500 ecat63InitMatlist(&mlist);
01501 ret=ecat63ReadMatlist(fp, &mlist);
01502 if(ret) {fclose(fp); return STATUS_NOMATLIST;}
01503 if(mlist.matrixNr<=0) {
01504 fclose(fp); ecat63EmptyMatlist(&mlist); return STATUS_INVALIDMATLIST;}
01505
01506 ecat63GatherMatlist(&mlist, 1, 1, 1, 1);
01507 ecat63SortMatlistByFrame(&mlist);
01508
01509
01510 ret=ecat63GetMatrixBlockSize(&mlist, &blkNr);
01511 if(ret) {ecat63EmptyMatlist(&mlist); fclose(fp); return ret;}
01512
01513 if(IMG_TEST>6) printf("allocating memory for %d blocks\n", blkNr);
01514 mdata=(char*)malloc(blkNr*MatBLKSIZE);
01515 if(mdata==NULL) {fclose(fp); ecat63EmptyMatlist(&mlist); return STATUS_NOMEMORY;}
01516
01517
01518 ret=0; seqplane=-1;
01519 for(m=0; m<mlist.matrixNr; m++) if(mlist.matdir[m].matstat==1) {
01520
01521 mat_numdoc(mlist.matdir[m].matnum, &matval);
01522 plane=matval.plane;
01523 if(main_header.num_frames>=main_header.num_gates) frame=matval.frame;
01524 else frame=matval.gate;
01525 if(frame!=frame_to_read) continue;
01526 seqplane++;
01527 if(img->planeNumber[seqplane]<1) {
01528 img->planeNumber[seqplane]=plane;
01529 } else if(img->planeNumber[seqplane]!=plane) {
01530 fclose(fp); ecat63EmptyMatlist(&mlist); free(mdata);
01531 return STATUS_MISSINGMATRIX;
01532 }
01533
01534
01535 if(IMG_TEST>4) printf("reading subheader for matrix %d,%d\n", frame, plane);
01536 if(main_header.file_type==IMAGE_DATA) {
01537 ret=ecat63ReadImageheader(fp, mlist.matdir[m].strtblk, &image_header);
01538 scale=image_header.quant_scale;
01539 if(image_header.ecat_calibration_fctr>0.0
01540 && fabs(main_header.calibration_factor/image_header.ecat_calibration_fctr-1.0)>0.0001)
01541 scale*=image_header.ecat_calibration_fctr;
01542 local_data_type=image_header.data_type;
01543 img->start[frame_index]=image_header.frame_start_time/1000.;
01544 img->end[frame_index]=img->start[frame_index]+image_header.frame_duration/1000.;
01545 img->mid[frame_index]=0.5*(img->start[frame_index]+img->end[frame_index]);
01546 if(image_header.decay_corr_fctr>1.0)
01547 img->decayCorrFactor[frame_index]=image_header.decay_corr_fctr;
01548 } else if(main_header.file_type==RAW_DATA) {
01549 ret=ecat63ReadScanheader(fp, mlist.matdir[m].strtblk, &scan_header);
01550 scale=scan_header.scale_factor;
01551 if(scan_header.loss_correction_fctr>0.0) scale*=scan_header.loss_correction_fctr;
01552 local_data_type=scan_header.data_type;
01553 img->start[frame_index]=scan_header.frame_start_time/1000.;
01554 img->end[frame_index]=img->start[frame_index]+scan_header.frame_duration/1000.;
01555 img->mid[frame_index]=0.5*(img->start[frame_index]+img->end[frame_index]);
01556 img->sampleDistance=10.0*scan_header.sample_distance;
01557 img->prompts[frame_index]=(float)scan_header.prompts;
01558 img->randoms[frame_index]=(float)scan_header.delayed;
01559 } else if(main_header.file_type==ATTN_DATA) {
01560 ret=ecat63ReadAttnheader(fp, mlist.matdir[m].strtblk, &attn_header);
01561 scale=attn_header.scale_factor;
01562 local_data_type=attn_header.data_type;
01563 } else if(main_header.file_type==NORM_DATA) {
01564 ret=ecat63ReadNormheader(fp, mlist.matdir[m].strtblk, &norm_header);
01565 scale=norm_header.scale_factor;
01566 local_data_type=norm_header.data_type;
01567 } else
01568 img->_dataType=-1;
01569 if(ret) {
01570 ecat63EmptyMatlist(&mlist); free(mdata); fclose(fp);
01571 return STATUS_NOSUBHEADER;
01572 }
01573 if(main_header.calibration_factor>0.0) scale*=main_header.calibration_factor;
01574 if(IMG_TEST>6) printf("scale=%e datatype=%d\n", scale, img->_dataType);
01575
01576 if(IMG_TEST>4) printf("reading matrix data\n");
01577 ret=ecat63ReadMatdata(fp, mlist.matdir[m].strtblk+1,
01578 mlist.matdir[m].endblk-mlist.matdir[m].strtblk,
01579 mdata, local_data_type);
01580 if(ret) {
01581 ecat63EmptyMatlist(&mlist); free(mdata); fclose(fp);
01582 return STATUS_MISSINGMATRIX;
01583 }
01584 if(IMG_TEST>4) printf("converting matrix data to floats\n");
01585 if(local_data_type==BYTE_TYPE) {
01586 for(yi=0, mptr=mdata; yi<img->dimy; yi++)
01587 for(xi=0; xi<img->dimx; xi++, mptr++) {
01588 img->m[seqplane][yi][xi][frame_index]=scale*(float)(*mptr);
01589 }
01590 } else if(local_data_type==VAX_I2 || local_data_type==SUN_I2) {
01591 for(yi=0, mptr=mdata; yi<img->dimy; yi++)
01592 for(xi=0; xi<img->dimx; xi++, mptr+=2) {
01593 sptr=(short int*)mptr;
01594 img->m[seqplane][yi][xi][frame_index]=scale*(float)(*sptr);
01595 }
01596 } else if(local_data_type==VAX_I4 || local_data_type==SUN_I4) {
01597 for(yi=0, mptr=mdata; yi<img->dimy; yi++)
01598 for(xi=0; xi<img->dimx; xi++, mptr+=4) {
01599 iptr=(int*)mptr;
01600 img->m[seqplane][yi][xi][frame_index]=1.0;
01601 }
01602 } else if(local_data_type==VAX_R4 || local_data_type==IEEE_R4) {
01603 for(yi=0, mptr=mdata; yi<img->dimy; yi++)
01604 for(xi=0; xi<img->dimx; xi++, mptr+=4) {
01605 memcpy(&img->m[seqplane][yi][xi][frame_index], mptr, 4);
01606 img->m[seqplane][yi][xi][frame_index]*=scale;
01607 }
01608 }
01609 }
01610 if(IMG_TEST>3) printf("end of matrices.\n");
01611
01612 free(mdata);
01613 ecat63EmptyMatlist(&mlist);
01614 fclose(fp);
01615
01616
01617 if(IMG_TEST>4) printf("last_seqplane := %d.\n", seqplane);
01618 if(seqplane<0) {imgSetStatus(img, STATUS_NOMATRIX); return STATUS_NOMATRIX;}
01619
01620
01621 if(seqplane+1 != img->dimz) {
01622 imgSetStatus(img, STATUS_MISSINGMATRIX);
01623 return STATUS_MISSINGMATRIX;
01624 }
01625
01626
01627 if(img->_dataType==VAX_I4 || img->_dataType==VAX_R4) img->_dataType=VAX_I2;
01628
01629
01630 imgSetStatus(img, STATUS_OK);
01631 return STATUS_OK;
01632 }
01633
01634
01635
01656 int imgWriteEcat63Frame(const char *fname, int frame_to_write, IMG *img, int frame_index) {
01657 IMG test_img;
01658 int ret=0, pxlNr, zi, xi, yi, matrixId;
01659 ECAT63_mainheader main_header;
01660 ECAT63_imageheader image_header;
01661 ECAT63_scanheader scan_header;
01662 void *sub_header=NULL;
01663 FILE *fp;
01664 float *fdata=NULL, *fptr;
01665
01666
01667 if(IMG_TEST) printf("\nimgWriteEcat63Frame(%s, %d, *img, %d)\n",
01668 fname, frame_to_write, frame_index);
01669 if(IMG_TEST>1) ECAT63_TEST=IMG_TEST-1;
01670
01671
01672
01673
01674 if(fname==NULL) return STATUS_FAULT;
01675 if(img==NULL) return STATUS_FAULT;
01676 if(img->status!=IMG_STATUS_OCCUPIED) return STATUS_FAULT;
01677 if(frame_to_write<0) return STATUS_FAULT;
01678 if(frame_index<0 || frame_index>=img->dimt) return STATUS_FAULT;
01679 if(img->_fileFormat!=IMG_E63) return STATUS_FAULT;
01680
01681
01682 memset(&main_header, 0, sizeof(ECAT63_mainheader));
01683 memset(&image_header,0, sizeof(ECAT63_imageheader));
01684 memset(&scan_header, 0, sizeof(ECAT63_scanheader));
01685 imgInit(&test_img);
01686
01687
01688
01689
01690
01691
01692
01693
01694 if(access(fname, 0) == -1) {
01695
01696
01697 imgSetEcat63MHeader(img, &main_header);
01698 if(IMG_TEST>2) {
01699 ecat63PrintMainheader(&main_header, stdout);
01700 }
01701 if(frame_to_write==0) frame_to_write=1;
01702 main_header.num_frames=frame_to_write;
01703
01704
01705 fp=ecat63Create(fname, &main_header); if(fp==NULL) return STATUS_NOWRITEPERM;
01706
01707 } else {
01708
01709
01710 ret=imgReadEcat63Header(fname, &test_img); if(ret!=0) return ret;
01711
01712 if(img->_fileFormat!=test_img._fileFormat || img->type!=test_img.type)
01713 return STATUS_WRONGFILETYPE;
01714
01715 if(img->dimz!=test_img.dimz || img->dimx!=test_img.dimx ||
01716 img->dimy!=test_img.dimy)
01717 return STATUS_VARMATSIZE;
01718 imgEmpty(&test_img);
01719
01720
01721 if((fp=fopen(fname, "r+b"))==NULL) return STATUS_NOWRITEPERM;
01722
01723
01724 if((ret=ecat63ReadMainheader(fp, &main_header))!=0) return STATUS_NOMAINHEADER;
01725 if(frame_to_write==0) frame_to_write=main_header.num_frames+1;
01726 if(main_header.num_frames<frame_to_write)
01727 main_header.num_frames=frame_to_write;
01728 if((ret=ecat63WriteMainheader(fp, &main_header))!=0) return STATUS_NOWRITEPERM;
01729
01730 }
01731 if(IMG_TEST>2) {
01732 printf("frame_to_write := %d\n", frame_to_write);
01733 }
01734
01735
01736 pxlNr=img->dimx*img->dimy*img->dimz;
01737 fdata=(float*)malloc(pxlNr*sizeof(float));
01738 if(fdata==NULL) {fclose(fp); return STATUS_NOMEMORY;}
01739
01740
01741 if(img->type==IMG_TYPE_RAW) sub_header=(void*)&scan_header;
01742 else if(img->type==IMG_TYPE_IMAGE) sub_header=&image_header;
01743 else {fclose(fp); return STATUS_FAULT;}
01744 imgSetEcat63SHeader(img, sub_header);
01745
01746
01747 fptr=fdata;
01748 for(zi=0; zi<img->dimz; zi++)
01749 for(yi=0; yi<img->dimy; yi++) for(xi=0; xi<img->dimx; xi++)
01750 *fptr++=img->m[zi][yi][xi][frame_index];
01751
01752
01753 fptr=fdata; ret=-1;
01754 for(zi=0; zi<img->dimz; zi++, fptr+=img->dimx*img->dimy) {
01755
01756 matrixId=mat_numcod(frame_to_write, img->planeNumber[zi], 1, 0, 0);
01757 if(img->type==IMG_TYPE_RAW) {
01758 scan_header.frame_start_time=
01759 (int)temp_roundf(1000.*img->start[frame_index]);
01760 scan_header.frame_duration=
01761 (int)temp_roundf(1000.*(img->end[frame_index]-img->start[frame_index]));
01762 scan_header.prompts=temp_roundf(img->prompts[frame_index]);
01763 scan_header.delayed=temp_roundf(img->randoms[frame_index]);
01764 ret=ecat63WriteScanMatrix(fp, matrixId, &scan_header, fptr);
01765 } else {
01766 image_header.frame_start_time=
01767 (int)temp_roundf(1000.*img->start[frame_index]);
01768 image_header.frame_duration=
01769 (int)temp_roundf(1000.*(img->end[frame_index]-img->start[frame_index]));
01770 if(img->decayCorrected)
01771 image_header.decay_corr_fctr=img->decayCorrFactor[frame_index];
01772 ret=ecat63WriteImageMatrix(fp, matrixId, &image_header, fptr);
01773 }
01774 if(ret!=0) break;
01775 }
01776 free(fdata); fclose(fp);
01777 if(ret) return STATUS_DISKFULL;
01778
01779 return STATUS_OK;
01780 }
01781
01782
01783
01790 void imgSetEcat63SHeader(IMG *img, void *h) {
01791 ECAT63_imageheader *image_header;
01792 ECAT63_scanheader *scan_header;
01793
01794 if(img->type==IMG_TYPE_RAW) {
01795 scan_header=(ECAT63_scanheader*)h;
01796 scan_header->data_type=VAX_I2;
01797 scan_header->dimension_1=img->dimx;
01798 scan_header->dimension_2=img->dimy;
01799 scan_header->frame_duration_sec=1;
01800 scan_header->scale_factor=1.0;
01801 scan_header->frame_start_time=0;
01802 scan_header->frame_duration=1000;
01803 scan_header->loss_correction_fctr=1.0;
01804 } else {
01805 image_header=(ECAT63_imageheader*)h;
01806 image_header->data_type=VAX_I2;
01807 image_header->num_dimensions=2;
01808 image_header->dimension_1=img->dimx;
01809 image_header->dimension_2=img->dimy;
01810 image_header->recon_scale=img->zoom;
01811 image_header->quant_scale=1.0;
01812 image_header->slice_width=img->sizez/10.;
01813 image_header->pixel_size=img->sizex/10.;
01814 image_header->frame_start_time=0;
01815 image_header->frame_duration=1000;
01816 image_header->plane_eff_corr_fctr=1.0;
01817 image_header->decay_corr_fctr=1.0;
01818 image_header->loss_corr_fctr=1.0;
01819 image_header->ecat_calibration_fctr=1.0;
01820 image_header->well_counter_cal_fctr=1.0;
01821 image_header->quant_units=imgUnitToEcat6(img);
01822 }
01823 }
01824
01825
01826