vol.c
Go to the documentation of this file.
00001 /******************************************************************************
00002   Copyright (c) 2003,2004 by Turku PET Centre
00003 
00004   vol.c
00005   
00006   Procedures for
00007   - storing and processing 3D PET image volume data (no time information).
00008   Depends on IMG functions in library.
00009   
00010   Version:
00011   2003-12-18 Vesa Oikonen
00012     Added 3D structures VOL and SVOL and related functions.
00013   2004-01-29 VO
00014     Added functions vol2img() and svol2img.
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   /*  0 */ "ok",
00034   /*  1 */ "fault in calling routine",
00035   /*  2 */ "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   /* Free up memory */
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   /* Set variables */
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   /* Set status */
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   /* Free up memory */
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   /* Set variables */
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   /* Set status */
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   /* Check arguments */
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   /* Check if correct volume size is already allocated */
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); /* that's it */
00173     } else {
00174       volEmpty(vol);
00175     }
00176   }
00177   /* Allocate memory for volume data */
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   /* Set data pointers */
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   /* Ok */
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   /* Check arguments */
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   /* Check if correct volume size is already allocated */
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); /* that's it */
00238     } else {
00239       svolEmpty(svol);
00240     }
00241   }
00242   /* Allocate memory for volume data */
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   /* Set data pointers */
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   /* Ok */
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   /* Check input */
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   /* Allocate memory (if needed) for volume */
00297   ret=volAllocate(vol, img->dimz, img->dimy, img->dimx);
00298   if(ret) return(ret);
00299 
00300   /* Copy data */
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   /* Check input */
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   /* Allocate memory (if needed) for volume */
00336   ret=svolAllocate(svol, img->dimz, img->dimy, img->dimx);
00337   if(ret) return(ret);
00338 
00339   /* Copy data */
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   /* Check input */
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   /* Copy data */
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   /* Check input */
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   /* Copy data */
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 /*****************************************************************************/