00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include <stdio.h>
00018 #include <stdlib.h>
00019 #include <math.h>
00020 #include <string.h>
00021 #include <time.h>
00022
00023 #include "petc99.h"
00024
00025 #include "include/img.h"
00026 #include "include/imgmax.h"
00027 #include "include/vol.h"
00028
00032 char *_volStatusMessage[] = {
00033 "ok",
00034 "fault in calling routine",
00035 "out of memory"
00036 };
00037
00038
00039
00045 void volInit(VOL *vol) {
00046 if(VOL_TEST) printf("volInit()\n");
00047 memset(vol, 0, sizeof(VOL));
00048 vol->status=IMG_STATUS_INITIALIZED;
00049 vol->statmsg=_volStatusMessage[0];
00050 vol->orientation=0;
00051 vol->dimx=vol->dimy=vol->dimz=0;
00052 vol->sizex=vol->sizey=vol->sizez=0;
00053 vol->v=(float***)NULL;
00054 vol->voxel=(float*)NULL;
00055 vol->column=(float*)NULL;
00056 vol->row=(float**)NULL;
00057 vol->plane=(float***)NULL;
00058 }
00059
00060
00061
00068 void svolInit(SVOL *svol) {
00069 if(VOL_TEST) printf("svolInit()\n");
00070 memset(svol, 0, sizeof(SVOL));
00071 svol->status=IMG_STATUS_INITIALIZED;
00072 svol->statmsg=_volStatusMessage[0];
00073 svol->orientation=0;
00074 svol->dimx=svol->dimy=svol->dimz=0;
00075 svol->sizex=svol->sizey=svol->sizez=0;
00076 svol->scale_factor=1.0;
00077 svol->v=(short int***)NULL;
00078 svol->voxel=(short int*)NULL;
00079 svol->column=(short int*)NULL;
00080 svol->row=(short int**)NULL;
00081 svol->plane=(short int***)NULL;
00082 }
00083
00084
00085
00091 void volEmpty(VOL *vol) {
00092 if(VOL_TEST) printf("volEmpty()\n");
00093 if(vol->status<IMG_STATUS_OCCUPIED) return;
00094
00095 if(vol->_vxl!=NULL) free(vol->_vxl);
00096 if(vol->_col!=NULL) free(vol->_col);
00097 if(vol->_row!=NULL) free(vol->_row);
00098 if(vol->_pln!=NULL) free(vol->_pln);
00099
00100 vol->statmsg=_volStatusMessage[0];
00101 vol->orientation=0;
00102 vol->dimx=vol->dimy=vol->dimz=0;
00103 vol->sizex=vol->sizey=vol->sizez=0;
00104 vol->v=(float***)NULL;
00105 vol->voxel=(float*)NULL;
00106 vol->column=(float*)NULL;
00107 vol->row=(float**)NULL;
00108 vol->plane=(float***)NULL;
00109
00110 vol->status=IMG_STATUS_INITIALIZED;
00111 }
00112
00113
00114
00120 void svolEmpty(SVOL *svol) {
00121 if(VOL_TEST) printf("svolEmpty()\n");
00122 if(svol->status<IMG_STATUS_OCCUPIED) return;
00123
00124 if(svol->_vxl!=NULL) free(svol->_vxl);
00125 if(svol->_col!=NULL) free(svol->_col);
00126 if(svol->_row!=NULL) free(svol->_row);
00127 if(svol->_pln!=NULL) free(svol->_pln);
00128
00129 svol->statmsg=_volStatusMessage[0];
00130 svol->orientation=0;
00131 svol->dimx=svol->dimy=svol->dimz=0;
00132 svol->sizex=svol->sizey=svol->sizez=0;
00133 svol->scale_factor=1.0;
00134 svol->v=(short int***)NULL;
00135 svol->voxel=(short int*)NULL;
00136 svol->column=(short int*)NULL;
00137 svol->row=(short int**)NULL;
00138 svol->plane=(short int***)NULL;
00139
00140 svol->status=IMG_STATUS_INITIALIZED;
00141 }
00142
00143
00144
00156 int volAllocate(VOL *vol, int planes, int rows, int columns) {
00157 unsigned short int zi, ri;
00158 int vxlNr, vi;
00159 float **rptr, *cptr;
00160
00161 if(VOL_TEST) printf("voiAllocate(*vol, %d, %d, %d)\n", planes, rows, columns);
00162
00163 if(vol==NULL) return(1); else vol->statmsg=_volStatusMessage[1];
00164 if(vol->status==IMG_STATUS_UNINITIALIZED) return(1);
00165 if(planes<1 || rows<1 || columns<1) return(2);
00166 vxlNr=planes*rows*columns;
00167
00168
00169 if(vol->status>=IMG_STATUS_OCCUPIED) {
00170 if(planes==vol->dimz && rows==vol->dimy && columns==vol->dimx) {
00171 for(vi=0; vi<vxlNr; vi++) vol->_vxl[vi]=0;
00172 return(0);
00173 } else {
00174 volEmpty(vol);
00175 }
00176 }
00177
00178 vol->_pln=(float***)malloc(planes*sizeof(float**));
00179 if(vol->_pln==NULL) {
00180 return(5);}
00181 vol->_row=(float**)malloc(planes*rows*sizeof(float*));
00182 if(vol->_row==NULL) {
00183 free(vol->_pln); return(6);}
00184 vol->_col=vol->_vxl=(float*)calloc(planes*rows*columns, sizeof(float));
00185 if(vol->_vxl==NULL) {
00186 free(vol->_pln); free(vol->_row); return(8);
00187 }
00188
00189 rptr=vol->_row; cptr=vol->_col;
00190 for(zi=0; zi<planes; zi++) {
00191 vol->_pln[zi]=rptr;
00192 for(ri=0; ri<rows; ri++) {
00193 *rptr++=cptr; cptr+=columns;
00194 }
00195 }
00196 vol->v=vol->_pln;
00197 vol->plane=vol->_pln;
00198 vol->column=vol->_col;
00199 vol->row=vol->_row;
00200 vol->voxel=vol->_vxl;
00201
00202 vol->dimz=planes; vol->dimy=rows; vol->dimx=columns;
00203 vol->statmsg=_volStatusMessage[0];
00204 vol->status=IMG_STATUS_OCCUPIED;
00205 return(0);
00206 }
00207
00208
00209
00221 int svolAllocate(SVOL *svol, int planes, int rows, int columns) {
00222 unsigned short int zi, ri;
00223 int vxlNr, vi;
00224 short int **rptr, *cptr;
00225
00226 if(VOL_TEST) printf("svoiAllocate(*svol, %d, %d, %d)\n", planes, rows, columns);
00227
00228 if(svol==NULL) return(1); else svol->statmsg=_volStatusMessage[1];
00229 if(svol->status==IMG_STATUS_UNINITIALIZED) return(1);
00230 if(planes<1 || rows<1 || columns<1) return(2);
00231 vxlNr=planes*rows*columns;
00232
00233
00234 if(svol->status>=IMG_STATUS_OCCUPIED) {
00235 if(planes==svol->dimz && rows==svol->dimy && columns==svol->dimx) {
00236 for(vi=0; vi<vxlNr; vi++) svol->_vxl[vi]=0;
00237 return(0);
00238 } else {
00239 svolEmpty(svol);
00240 }
00241 }
00242
00243 svol->_pln=(short int***)malloc(planes*sizeof(short int**));
00244 if(svol->_pln==NULL) {
00245 return(5);}
00246 svol->_row=(short int**)malloc(planes*rows*sizeof(short int*));
00247 if(svol->_row==NULL) {
00248 free(svol->_pln); return(6);}
00249 svol->_col=svol->_vxl=(short int*)calloc(planes*rows*columns, sizeof(short int));
00250 if(svol->_vxl==NULL) {
00251 free(svol->_pln); free(svol->_row); return(8);
00252 }
00253
00254 rptr=svol->_row; cptr=svol->_col;
00255 for(zi=0; zi<planes; zi++) {
00256 svol->_pln[zi]=rptr;
00257 for(ri=0; ri<rows; ri++) {
00258 *rptr++=cptr; cptr+=columns;
00259 }
00260 }
00261 svol->v=svol->_pln;
00262 svol->plane=svol->_pln;
00263 svol->column=svol->_col;
00264 svol->row=svol->_row;
00265 svol->voxel=svol->_vxl;
00266
00267 svol->dimz=planes; svol->dimy=rows; svol->dimx=columns;
00268 svol->statmsg=_volStatusMessage[0];
00269 svol->status=IMG_STATUS_OCCUPIED;
00270 return(0);
00271 }
00272
00273
00274
00284 int img2vol(IMG *img, VOL *vol, int frame) {
00285 int ret;
00286 unsigned short int zi, yi, xi, fi;
00287
00288 if(VOL_TEST) printf("img2vol(img, %d, vol)\n", frame);
00289
00290 if(vol==NULL) return(1);
00291 vol->statmsg=_volStatusMessage[1];
00292 if(img->status!=IMG_STATUS_OCCUPIED) return(1);
00293 if(frame<1 || img->dimt<frame) return(2);
00294 if(vol->status==IMG_STATUS_UNINITIALIZED) return(1);
00295
00296
00297 ret=volAllocate(vol, img->dimz, img->dimy, img->dimx);
00298 if(ret) return(ret);
00299
00300
00301 fi=frame-1;
00302 vol->orientation=img->orientation;
00303 vol->sizex=img->sizex; vol->sizey=img->sizey; vol->sizez=img->sizez;
00304 for(zi=0; zi<vol->dimz; zi++)
00305 for(yi=0; yi<vol->dimy; yi++)
00306 for(xi=0; xi<vol->dimx; xi++)
00307 vol->v[zi][yi][xi]=img->m[zi][yi][xi][fi];
00308
00309 return(0);
00310 }
00311
00312
00313
00322 int img2svol(IMG *img, SVOL *svol, int frame) {
00323 int ret;
00324 unsigned short int zi, yi, xi, fi;
00325 float fmin, fmax, g;
00326
00327 if(VOL_TEST) printf("img2svol(img, %d, svol)\n", frame);
00328
00329 if(svol==NULL) return(1);
00330 svol->statmsg=_volStatusMessage[1];
00331 if(img->status!=IMG_STATUS_OCCUPIED) return(1);
00332 if(frame<1 || img->dimt<frame) return(2);
00333 if(svol->status==IMG_STATUS_UNINITIALIZED) return(1);
00334
00335
00336 ret=svolAllocate(svol, img->dimz, img->dimy, img->dimx);
00337 if(ret) return(ret);
00338
00339
00340 fi=frame-1;
00341 svol->orientation=img->orientation;
00342 svol->sizex=img->sizex; svol->sizey=img->sizey; svol->sizez=img->sizez;
00343 ret=imgFrameMinMax(img, frame, &fmin, &fmax); if(ret) return(10+ret);
00344 if(fabs(fmin)>fabs(fmax)) g=fabs(fmin); else g=fabs(fmax);
00345 if(g!=0) g=32766./g; else g=1.0;
00346 for(zi=0; zi<svol->dimz; zi++)
00347 for(yi=0; yi<svol->dimy; yi++)
00348 for(xi=0; xi<svol->dimx; xi++)
00349 svol->v[zi][yi][xi]=(short int)temp_roundf(g*img->m[zi][yi][xi][fi]);
00350 svol->scale_factor=1.0/g;
00351
00352 return(0);
00353 }
00354
00355
00356
00367 int vol2img(VOL *vol, IMG *img, int frame) {
00368 unsigned short int zi, yi, xi, fi;
00369
00370 if(VOL_TEST) printf("vol2img(vol, img, %d)\n", frame);
00371
00372 if(vol==NULL || vol->status!=IMG_STATUS_OCCUPIED) return(1);
00373 vol->statmsg=_volStatusMessage[1];
00374 if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) return(1);
00375 if(frame<1 || img->dimt<frame) return(2);
00376 if(img->dimx!=vol->dimx || img->dimy!=vol->dimy) return(3);
00377 if(img->dimz!=vol->dimz) return(4);
00378
00379
00380 fi=frame-1;
00381 for(zi=0; zi<vol->dimz; zi++)
00382 for(yi=0; yi<vol->dimy; yi++)
00383 for(xi=0; xi<vol->dimx; xi++)
00384 img->m[zi][yi][xi][fi]=vol->v[zi][yi][xi];
00385
00386 return(0);
00387 }
00388
00389
00390
00401 int svol2img(SVOL *svol, IMG *img, int frame) {
00402 unsigned short int zi, yi, xi, fi;
00403
00404 if(VOL_TEST) printf("svol2img(svol, img, %d)\n", frame);
00405
00406 if(svol==NULL || svol->status!=IMG_STATUS_OCCUPIED) return(1);
00407 svol->statmsg=_volStatusMessage[1];
00408 if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) return(1);
00409 if(frame<1 || img->dimt<frame) return(2);
00410 if(img->dimx!=svol->dimx || img->dimy!=svol->dimy) return(3);
00411 if(img->dimz!=svol->dimz) return(4);
00412
00413
00414 fi=frame-1;
00415 for(zi=0; zi<svol->dimz; zi++)
00416 for(yi=0; yi<svol->dimy; yi++)
00417 for(xi=0; xi<svol->dimx; xi++)
00418 img->m[zi][yi][xi][fi]=(svol->scale_factor)*(float)svol->v[zi][yi][xi];
00419
00420 return(0);
00421 }
00422
00423
00424
00431 void volInfo(VOL *vol, FILE *fp) {
00432 if(VOL_TEST) printf("volInfo()\n");
00433 if(vol->status<=IMG_STATUS_UNINITIALIZED) {
00434 fprintf(fp, "Volume data is not initialized.\n"); return;}
00435 if(vol->status==IMG_STATUS_INITIALIZED) {
00436 fprintf(fp, "Volume data is initialized but empty.\n"); return;}
00437 if(vol->status==IMG_STATUS_ERROR) fprintf(stdout, "Volume data has errors.\n");
00438 fprintf(fp, "Volume status: %s\n", vol->statmsg);
00439 fprintf(fp, "Patient orientation: %d\n", vol->orientation);
00440 fprintf(fp, "Voxel sizes (x, y, z): %g %g %g mm\n",
00441 vol->sizex, vol->sizey, vol->sizez);
00442 fprintf(fp, "Dimensions (x, y, z): %d %d %d\n",
00443 vol->dimx, vol->dimy, vol->dimz);
00444 return;
00445 }
00446
00447
00448
00455 void svolInfo(SVOL *svol, FILE *fp) {
00456 if(VOL_TEST) printf("svolInfo()\n");
00457 if(svol->status<=IMG_STATUS_UNINITIALIZED) {
00458 fprintf(fp, "Volume data is not initialized.\n"); return;}
00459 if(svol->status==IMG_STATUS_INITIALIZED) {
00460 fprintf(fp, "Volume data is initialized but empty.\n"); return;}
00461 if(svol->status==IMG_STATUS_ERROR) fprintf(stdout, "Volume data has errors.\n");
00462 fprintf(fp, "Volume status: %s\n", svol->statmsg);
00463 fprintf(fp, "Patient orientation: %d\n", svol->orientation);
00464 fprintf(fp, "Voxel sizes (x, y, z): %g %g %g mm\n",
00465 svol->sizex, svol->sizey, svol->sizez);
00466 fprintf(fp, "Dimensions (x, y, z): %d %d %d\n",
00467 svol->dimx, svol->dimy, svol->dimz);
00468 fprintf(fp, "Scale factor: %g\n", svol->scale_factor);
00469 return;
00470 }
00471
00472
00473
00480 void volContents(VOL *vol, VOL_RANGE r, FILE *fp) {
00481 int zi, yi, xi;
00482
00483 if(vol->status!=IMG_STATUS_OCCUPIED) return;
00484 if(r.z1<1 || r.y1<1 || r.x1<1) return;
00485 if(r.z2<r.z1 || r.y2<r.y1 || r.x2<r.x1) return;
00486 if(r.z2>vol->dimz || r.y2>vol->dimy || r.x2>vol->dimx) return;
00487
00488 for(zi=r.z1-1; zi<r.z2; zi++) {
00489 fprintf(fp, "pl=%03d ", zi+1);
00490 for(xi=r.x1-1; xi<r.x2; xi++) fprintf(fp, " x=%05d", xi+1);
00491 fprintf(fp, "\n");
00492 for(yi=r.y1-1; yi<r.y2; yi++) {
00493 fprintf(fp, "y=%05d", yi+1);
00494 for(xi=r.x1-1; xi<r.x2; xi++)
00495 fprintf(fp, " %7.3f", vol->v[zi][yi][xi]);
00496 fprintf(fp, "\n");
00497 }
00498 }
00499 }
00500
00501
00502
00513 int volMax(VOL *vol, VOL_RANGE r, VOL_PIXEL *p, float *maxv) {
00514 int zi, yi, xi;
00515
00516 if(vol->status!=IMG_STATUS_OCCUPIED) return(1);
00517 if(r.z1<1 || r.y1<1 || r.x1<1) return(2);
00518 if(r.z2<r.z1 || r.y2<r.y1 || r.x2<r.x1) return(3);
00519 if(r.z2>vol->dimz || r.y2>vol->dimy || r.x2>vol->dimx) return(4);
00520
00521 zi=r.z1-1; yi=r.y1-1; xi=r.x1-1; *maxv=vol->v[zi][yi][xi];
00522 if(p!=NULL) {p->z=zi+1; p->y=yi+1; p->x=xi+1;}
00523 for(zi=r.z1-1; zi<r.z2; zi++) {
00524 for(yi=r.y1-1; yi<r.y2; yi++) {
00525 for(xi=r.x1-1; xi<r.x2; xi++) if(*maxv<vol->v[zi][yi][xi]) {
00526 *maxv=vol->v[zi][yi][xi];
00527 if(p!=NULL) {p->z=zi+1; p->y=yi+1; p->x=xi+1;}
00528 }
00529 }
00530 }
00531 return(0);
00532 }
00533
00534
00535
00545 int volAvg(VOL *vol, VOL_RANGE r, float *avg) {
00546 int zi, yi, xi, n=0;
00547
00548 if(vol->status!=IMG_STATUS_OCCUPIED) return(1);
00549 if(r.z1<1 || r.y1<1 || r.x1<1) return(2);
00550 if(r.z2<r.z1 || r.y2<r.y1 || r.x2<r.x1) return(3);
00551 if(r.z2>vol->dimz || r.y2>vol->dimy || r.x2>vol->dimx) return(4);
00552
00553 *avg=0.0;
00554 for(zi=r.z1-1; zi<r.z2; zi++) {
00555 for(yi=r.y1-1; yi<r.y2; yi++) {
00556 for(xi=r.x1-1; xi<r.x2; xi++) {
00557 *avg+=vol->v[zi][yi][xi]; n++;
00558 }
00559 }
00560 }
00561 if(n>0) *avg/=(float)n;
00562 return(0);
00563 }
00564
00565
00566