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
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100 #include <stdio.h>
00101 #include <stdlib.h>
00102 #include <math.h>
00103 #include <string.h>
00104 #include <time.h>
00105
00106 #include "studynr.h"
00107 #include "halflife.h"
00108
00109 #include "include/imgunit.h"
00110 #include "include/img.h"
00111
00113 static const char *imgmsg[] = {
00114 "ok",
00115 "fault in calling routine",
00116 "out of memory",
00117 "cannot open file",
00118 "unknown file format",
00119 "unsupported file type",
00120 "missing matrix/matrices",
00121 "no permission to write",
00122 "disk full?",
00123 "cannot read matrix list",
00124 "invalid matrix list",
00125 "variable matrix size",
00126 "cannot read mainheader",
00127 "cannot read subheader",
00128 "cannot read matrix",
00129 "axial compression is not supported",
00130 "image datafile does not exist",
00131 "header file does not exist",
00132 "invalid header contents",
00133 "cannot read image data",
00134 "cannot read sif data",
00135 "wrong sif data",
00136 "cannot write image datafile",
00137 "cannot write header file",
00138 "wrong file type",
00139 "cannot erase file",
00140 "cannot read data",
00141 "cannot write data",
00142 "polar map is not supported",
00143 "invalid polar map",
00144 0
00145 };
00146 #if 0
00147
00148 static const char *_imgStatusMessage[] =
00149 {
00150 "ok",
00151 "fault in calling routine",
00152 "out of memory"
00153 };
00154 #endif
00155
00156
00157
00163 void imgInit(IMG *image) {
00164 int i = 0;
00165 if(IMG_TEST) printf("imgInit()\n");
00166 memset(image, 0, sizeof(IMG));
00167
00168 image->status=IMG_STATUS_INITIALIZED;
00169 imgSetStatus(image, STATUS_OK);
00170 image->type=0;
00171 image->unit=0;
00172 image->zoom=0.0;
00173 image->radiopharmaceutical[0]=(char)0;
00174 image->isotopeHalflife=0.0;
00175 image->decayCorrected=(char)0;
00176 image->unit=0;
00177 image->scanStart=0;
00178 image->orientation=0;
00179 image->axialFOV=image->transaxialFOV=image->sampleDistance=0.0;
00180 image->studyNr[0]=image->patientName[0]=(char)0;
00181 image->sizex=image->sizey=image->sizez=0;
00182 image->_dataType=0;
00183 image->_fileFormat=0;
00184 image->scanner=0;
00185 image->polarmap_num_rings=0;
00186 for(i=0; i<MAX_POLARMAP_NUM_RINGS; i++) {
00187 image->polarmap_sectors_per_ring[i]=0;
00188 image->polarmap_ring_position[i]=0.0;
00189 image->polarmap_ring_angle[i]=0;
00190 }
00191 image->polarmap_start_angle=0;
00192 image->dimt=image->dimx=image->dimy=image->dimz=0;
00193 image->gapx=image->gapy=image->gapz=0.0;
00194 image->resolutionx=image->resolutiony=image->resolutionz=0.0;
00195 image->m=(float****)NULL;
00196 image->_header=(float*)NULL;
00197 image->pixel=(float*)NULL;
00198 image->column=(float**)NULL;
00199 image->row=(float***)NULL;
00200 image->plane=(float****)NULL;
00201 image->planeNumber=(int*)NULL;
00202 image->start=image->end=image->mid=(float*)NULL;
00203 image->isWeight=0;
00204 image->weight=image->sd=image->prompts=image->randoms=(float*)NULL;
00205 image->decayCorrFactor=(float*)NULL;
00206 image->errstatus=STATUS_OK;
00207 }
00208
00209
00210
00216 void imgEmpty(IMG *image) {
00217 int i = 0;
00218 if(IMG_TEST) printf("imgEmpty()\n");
00219 if(image->status<IMG_STATUS_OCCUPIED) return;
00220
00221 if(image->_pxl!=NULL) free(image->_pxl);
00222 if(image->_col!=NULL) free(image->_col);
00223 if(image->_row!=NULL) free(image->_row);
00224 if(image->_pln!=NULL) free(image->_pln);
00225 if(image->dimz>0) free(image->planeNumber);
00226 if(image->dimt>0) free(image->_header);
00227
00228 imgSetStatus(image, STATUS_OK);
00229 image->type=0;
00230 image->unit=0;
00231 image->zoom=0.0;
00232 image->radiopharmaceutical[0]=(char)0;
00233 image->isotopeHalflife=0.0;
00234 image->decayCorrected=(char)0;
00235 image->unit=0;
00236 image->scanStart=0;
00237 image->orientation=0;
00238 image->axialFOV=image->transaxialFOV=image->sampleDistance=0.0;
00239 image->studyNr[0]=image->patientName[0]=image->patientID[0]=(char)0;
00240 image->userProcessCode[0]=image->studyDescription[0]=(char)0;
00241 image->sizex=image->sizey=image->sizez=0;
00242 image->gapx=image->gapy=image->gapz=0.0;
00243 image->resolutionx=image->resolutiony=image->resolutionz=0.0;
00244 image->_dataType=0;
00245 image->_fileFormat=0;
00246 image->scanner=0;
00247 image->polarmap_num_rings=0;
00248 for(i=0; i<MAX_POLARMAP_NUM_RINGS; i++) {
00249 image->polarmap_sectors_per_ring[i]=0;
00250 image->polarmap_ring_position[i]=0.0;
00251 image->polarmap_ring_angle[i]=0;
00252 }
00253 image->polarmap_start_angle=0;
00254 image->dimt=image->dimx=image->dimy=image->dimz=0;
00255 image->m=(float****)NULL;
00256 image->_header=(float*)NULL;
00257 image->pixel=(float*)NULL;
00258 image->column=(float**)NULL;
00259 image->row=(float***)NULL;
00260 image->plane=(float****)NULL;
00261 image->planeNumber=(int*)NULL;
00262 image->start=image->end=image->mid=(float*)NULL;
00263 image->isWeight=0;
00264 image->weight=image->sd=(float*)NULL;
00265 image->decayCorrFactor=(float*)NULL;
00266
00267 image->status=IMG_STATUS_INITIALIZED;
00268 image->errstatus=STATUS_OK;
00269 }
00270
00271
00272
00285 int imgAllocate(IMG *image, int planes, int rows, int columns, int frames) {
00286 unsigned short int zi, ri, ci;
00287 float ***rptr, **cptr, *pptr;
00288
00289 if(IMG_TEST) printf("imgAllocate(*image, %d, %d, %d, %d)\n", planes, rows, columns, frames);
00290
00291 imgSetStatus(image, STATUS_FAULT);
00292 if(image->status==IMG_STATUS_UNINITIALIZED) return(1);
00293 if(planes<1 || rows<1 || columns<1 || frames<1) return(2);
00294 if(image->status>=IMG_STATUS_OCCUPIED) imgEmpty(image);
00295
00296 imgSetStatus(image, STATUS_NOMEMORY);
00297 image->_header=(float*)calloc(8*frames, sizeof(float));
00298 if(image->_header==NULL)
00299 return(3);
00300 image->planeNumber=(int*)calloc(planes, sizeof(int));
00301 if(image->planeNumber==NULL) {
00302 free(image->_header); return(4);}
00303
00304 image->_pln=(float****)malloc(planes*sizeof(float***));
00305 if(image->_pln==NULL) {
00306 free(image->_header); free(image->planeNumber); return(5);}
00307 image->_row=(float***)malloc(planes*rows*sizeof(float**));
00308 if(image->_row==NULL) {
00309 free(image->_header); free(image->planeNumber); free(image->_pln); return(6);}
00310 image->_col=(float**)malloc(planes*rows*columns*sizeof(float*));
00311 if(image->_col==NULL) {
00312 free(image->_header); free(image->planeNumber);
00313 free(image->_pln); free(image->_row); return(7);
00314 }
00315 image->_pxl=(float*)calloc(planes*rows*columns*frames, sizeof(float));
00316 if(image->_pxl==NULL) {
00317 free(image->_header); free(image->planeNumber);
00318 free(image->_pln); free(image->_row); free(image->_col); return(8);
00319 }
00320
00321 rptr=image->_row; cptr=image->_col; pptr=image->_pxl;
00322 for(zi=0; zi<planes; zi++) {
00323 image->_pln[zi]=rptr;
00324 for(ri=0; ri<rows; ri++) {
00325 *rptr++=cptr;
00326 for(ci=0; ci<columns; ci++) {
00327 *cptr++=pptr; pptr+=frames;
00328 }
00329 }
00330 }
00331 image->m=image->_pln;
00332 image->plane=image->_pln;
00333 image->column=image->_col;
00334 image->row=image->_row;
00335 image->pixel=image->_pxl;
00336
00337 image->start= image->_header+0*frames;
00338 image->end= image->_header+1*frames;
00339 image->mid= image->_header+2*frames;
00340 image->weight= image->_header+3*frames;
00341 image->sd= image->_header+4*frames;
00342 image->prompts= image->_header+5*frames;
00343 image->randoms= image->_header+6*frames;
00344 image->decayCorrFactor=image->_header+7*frames;
00345
00346 image->dimz=planes; image->dimy=rows; image->dimx=columns; image->dimt=frames;
00347 imgSetStatus(image, STATUS_OK);
00348 image->status=IMG_STATUS_OCCUPIED;
00349 return(0);
00350 }
00351
00352
00353
00357 int imgAllocateWithHeader(
00359 IMG *image,
00361 int planes,
00363 int rows,
00365 int columns,
00367 int frames,
00369 IMG *image_from
00370 ) {
00371 int ret;
00372 ret=imgAllocate(image, planes, rows, columns, frames); if(ret) return ret;
00373 ret=imgCopyhdr(image_from, image); return ret;
00374 }
00375
00376
00377
00384 char *imgStatus(int status_index) {
00385 int n=0;
00386 while(imgmsg[n]!=0) n++;
00387 if(status_index<0 || status_index>n-1) return((char*)imgmsg[STATUS_FAULT]);
00388 else return((char*)imgmsg[status_index]);
00389 }
00390
00391
00392
00399 void imgSetStatus(IMG *img, int status_index) {
00400 int n=0;
00401 while(imgmsg[n]!=0) n++;
00402 if(status_index<0 || status_index>n-1) img->errstatus=STATUS_FAULT;
00403 else img->errstatus=status_index;
00404 img->statmsg=imgmsg[img->errstatus];
00405 }
00406
00407
00408
00414 void imgInfo(IMG *image) {
00415 int i;
00416 char buf[64];
00417 struct tm *st;
00418
00419 if(IMG_TEST) printf("imgInfo()\n");
00420 if(image->status<=IMG_STATUS_UNINITIALIZED) {
00421 fprintf(stdout, "image_status := not initialized\n"); return;
00422 } else if(image->status==IMG_STATUS_INITIALIZED) {
00423 fprintf(stdout, "image_status := initialized but empty\n");
00424 } else if(image->status==IMG_STATUS_ERROR) {
00425 fprintf(stdout, "image_status := error\n");
00426 }
00427 fprintf(stdout, "image_error_status := %s\n", image->statmsg);
00428 fprintf(stdout, "image_type := %d\n", image->type);
00429 fprintf(stdout, "saved_data_type := %d\n", image->_dataType);
00430 fprintf(stdout, "file_format := %d\n", image->_fileFormat);
00431 fprintf(stdout, "scanner := %d\n", image->scanner);
00432 fprintf(stdout, "identification_code := %.*s\n", MAX_STUDYNR_LEN, image->studyNr);
00433 fprintf(stdout, "data_unit := %s\n", imgUnit((int)image->unit));
00434 fprintf(stdout, "image_zoom := %g\n", image->zoom);
00435 fprintf(stdout, "radiopharmaceutical := %.32s\n", image->radiopharmaceutical);
00436 fprintf(stdout, "isotope_halflife := %e [sec]\n", image->isotopeHalflife);
00437 st=localtime(&image->scanStart);
00438 if(st!=NULL) strftime(buf, 64, "%Y-%m-%d %H:%M:%S", st); else strcpy(buf, "1900-01-01 00:00:00");
00439 fprintf(stdout, "scan_start_time := %s\n", buf);
00440 fprintf(stdout, "patient_orientation := %d\n", image->orientation);
00441 fprintf(stdout, "FOV_axial := %g [mm]\n", image->axialFOV);
00442 fprintf(stdout, "FOV_transaxial := %g [mm]\n", image->transaxialFOV);
00443 fprintf(stdout, "sample_distance := %g [mm]\n", image->sampleDistance);
00444 fprintf(stdout, "pixel_size_x := %g [mm]\n", image->sizex);
00445 fprintf(stdout, "pixel_size_y := %g [mm]\n", image->sizey);
00446 fprintf(stdout, "pixel_size_z := %g [mm]\n", image->sizez);
00447 fprintf(stdout, "dimension_x := %d\n", image->dimx);
00448 fprintf(stdout, "dimension_y := %d\n", image->dimy);
00449 fprintf(stdout, "dimension_z := %d\n", image->dimz);
00450 fprintf(stdout, "dimension_t := %d\n", image->dimt);
00451
00452 fprintf(stdout, "polarmap_num_rings := %d\n", image->polarmap_num_rings);
00453 if(image->polarmap_num_rings>0) {
00454 fprintf(stdout, "polarmap_sectors_per_ring :=");
00455 for(i=0; i<image->polarmap_num_rings; i++)
00456 fprintf(stdout, " %d", image->polarmap_sectors_per_ring[i]);
00457 fprintf(stdout, "\n");
00458 fprintf(stdout, "polarmap_ring_position :=");
00459 for(i=0; i<image->polarmap_num_rings; i++)
00460 fprintf(stdout, " %g", image->polarmap_ring_position[i]);
00461 fprintf(stdout, "\n");
00462 fprintf(stdout, "polarmap_ring_angle :=");
00463 for(i=0; i<image->polarmap_num_rings; i++)
00464 fprintf(stdout, " %d", image->polarmap_ring_angle[i]);
00465 fprintf(stdout, "\n");
00466 fprintf(stdout, "polarmap_start_angle := %d\n", image->polarmap_start_angle);
00467 }
00468
00469 if(image->status==IMG_STATUS_INITIALIZED) return;
00470
00471 fprintf(stdout, "actual_plane_numbers := %d", image->planeNumber[0]);
00472 for(i=1; i<image->dimz; i++) fprintf(stdout, " %d", image->planeNumber[i]);
00473 fprintf(stdout, "\n");
00474 fprintf(stdout, "Frame times (sec):\n");
00475 for(i=0; i<image->dimt; i++) fprintf(stdout, " %e %e %e\n",
00476 image->start[i], image->end[i], image->mid[i]);
00477 if(image->isWeight) fprintf(stdout, "Frames are weighted.\n");
00478 else fprintf(stdout, "Frames are not weighted.\n");
00479 if(image->decayCorrected==1) {
00480 fprintf(stdout, "Decay correction factors for each frame:\n");
00481 for(i=0; i<image->dimt; i++)
00482 fprintf(stdout, "%03i %e\n", i+1, image->decayCorrFactor[i]);
00483 } else
00484 fprintf(stdout, "Image is not decay corrected.\n");
00485 return;
00486 }
00487
00488
00489
00501 int imgCopyhdr(IMG *image1, IMG *image2) {
00502 int i;
00503
00504 if(IMG_TEST) printf("imgCopyhdr()\n");
00505
00506 if(image1==NULL || image2==NULL) return(1);
00507 if(image1==image2) return(2);
00508
00509 image2->type=image1->type;
00510 image2->unit=image1->unit;
00511 strcpy(image2->studyNr, image1->studyNr);
00512 strcpy(image2->patientName, image1->patientName);
00513 strcpy(image2->patientID, image1->patientID);
00514 strcpy(image2->userProcessCode, image1->userProcessCode);
00515 strcpy(image2->studyDescription, image1->studyDescription);
00516 image2->zoom=image1->zoom;
00517 strcpy(image2->radiopharmaceutical, image1->radiopharmaceutical);
00518 image2->isotopeHalflife=image1->isotopeHalflife;
00519 image2->decayCorrected=image1->decayCorrected;
00520 image2->scanStart=image1->scanStart;
00521 image2->axialFOV=image1->axialFOV;
00522 image2->transaxialFOV=image1->transaxialFOV;
00523 image2->sampleDistance=image1->sampleDistance;
00524 image2->sizex=image1->sizex;
00525 image2->sizey=image1->sizey;
00526 image2->sizez=image1->sizez;
00527 image2->gapx=image1->gapx;
00528 image2->gapy=image1->gapy;
00529 image2->gapz=image1->gapz;
00530 image2->resolutionx=image1->resolutionx;
00531 image2->resolutiony=image1->resolutiony;
00532 image2->resolutionz=image1->resolutionz;
00533 image2->_dataType=image1->_dataType;
00534 image2->_fileFormat=image1->_fileFormat;
00535 image2->orientation=image1->orientation;
00536 image2->scanner=image1->scanner;
00537 image2->polarmap_num_rings=image1->polarmap_num_rings;
00538 for(i=0; i<MAX_POLARMAP_NUM_RINGS; i++) {
00539 image2->polarmap_sectors_per_ring[i]=image1->polarmap_sectors_per_ring[i];
00540 image2->polarmap_ring_position[i]=image1->polarmap_ring_position[i];
00541 image2->polarmap_ring_angle[i]=image1->polarmap_ring_angle[i];
00542 }
00543 image2->polarmap_start_angle=image1->polarmap_start_angle;
00544 if(image1->dimz==image2->dimz) for(i=0; i<image1->dimz; i++)
00545 image2->planeNumber[i]=image1->planeNumber[i];
00546 if(image1->dimt==image2->dimt) for(i=0; i<image1->dimt; i++) {
00547 image2->start[i]=image1->start[i]; image2->end[i]=image1->end[i];
00548 image2->mid[i]=image1->mid[i];
00549 image2->weight[i]=image1->weight[i]; image2->sd[i]=image1->sd[i];
00550 image2->prompts[i]=image1->prompts[i];
00551 image2->randoms[i]=image1->randoms[i];
00552 image2->decayCorrFactor[i]=image1->decayCorrFactor[i];
00553 }
00554 image2->isWeight=image1->isWeight;
00555 return(0);
00556 }
00557
00558
00559
00570 int imgExtractRange(IMG *img1, IMG_RANGE r, IMG *img2) {
00571 int zi, yi, xi, fi, zj, yj, xj, fj;
00572
00573 if(IMG_TEST) {
00574 printf("imgExtractRange(*img1, r, *img2)\n");
00575 printf(" z=[%d,%d] y=[%d,%d] x=[%d,%d] f=[%d,%d]\n",
00576 r.z1, r.z2, r.y1, r.y2, r.x1, r.x2, r.f1, r.f2);
00577 }
00578
00579 if(img2==NULL) return(1); else imgSetStatus(img2, STATUS_FAULT);
00580 if(img1->status!=IMG_STATUS_OCCUPIED) return(1);
00581 if(img2->status==IMG_STATUS_UNINITIALIZED) return(1);
00582 if(r.z1<1 || r.z2>img1->dimz || r.z1>r.z2) return(1);
00583 if(r.y1<1 || r.y2>img1->dimy || r.y1>r.y2) return(1);
00584 if(r.x1<1 || r.x2>img1->dimx || r.x1>r.x2) return(1);
00585 if(r.f1<1 || r.f2>img1->dimt || r.f1>r.f2) return(1);
00586
00587
00588 imgSetStatus(img2, STATUS_NOMEMORY);
00589 zi=r.z2-r.z1+1; yi=r.y2-r.y1+1; xi=r.x2-r.x1+1; fi=r.f2-r.f1+1;
00590 if(img2->status>=IMG_STATUS_OCCUPIED)
00591 if(img2->dimz!=zi || img2->dimy!=yi || img2->dimx!=xi || img2->dimt!=fi)
00592 imgEmpty(img2);
00593 if(img2->status!=IMG_STATUS_OCCUPIED) {
00594 if(imgAllocate(img2, zi, yi, xi, fi)!=0) return(2);
00595 }
00596
00597
00598 imgCopyhdr(img1, img2);
00599 for(fi=r.f1-1, fj=0; fi<r.f2; fi++, fj++) {
00600 img2->start[fj]=img1->start[fi];
00601 img2->end[fj]=img1->end[fi];
00602 img2->mid[fj]=img1->mid[fi];
00603 img2->weight[fj]=img1->weight[fi];
00604 img2->sd[fj]=img1->sd[fi];
00605 img2->prompts[fj]=img1->prompts[fi];
00606 img2->randoms[fj]=img1->randoms[fi];
00607 img2->decayCorrFactor[fj]=img1->decayCorrFactor[fi];
00608 }
00609 for(zi=r.z1-1, zj=0; zi<r.z2; zi++, zj++)
00610 img2->planeNumber[zj]=img1->planeNumber[zi];
00611 for(zi=r.z1-1, zj=0; zi<r.z2; zi++, zj++)
00612 for(yi=r.y1-1, yj=0; yi<r.y2; yi++, yj++)
00613 for(xi=r.x1-1, xj=0; xi<r.x2; xi++, xj++)
00614 for(fi=r.f1-1, fj=0; fi<r.f2; fi++, fj++)
00615 img2->m[zj][yj][xj][fj]=img1->m[zi][yi][xi][fi];
00616
00617 imgSetStatus(img2, STATUS_OK);
00618 return(0);
00619 }
00620
00621
00622