vol.c
Go to the documentation of this file.
1 /******************************************************************************
2  Copyright (c) 2003,2004 by Turku PET Centre
3 
4  vol.c
5 
6  Procedures for
7  - storing and processing 3D PET image volume data (no time information).
8  Depends on IMG functions in library.
9 
10  Version:
11  2003-12-18 Vesa Oikonen
12  Added 3D structures VOL and SVOL and related functions.
13  2004-01-29 VO
14  Added functions vol2img() and svol2img.
15 
16 ******************************************************************************/
17 #include <stdio.h>
18 #include <stdlib.h>
19 #include <math.h>
20 #include <string.h>
21 #include <time.h>
22 /*****************************************************************************/
23 #include "petc99.h"
24 /*****************************************************************************/
25 #include "include/img.h"
26 #include "include/imgmax.h"
27 #include "include/vol.h"
28 /*****************************************************************************/
32 char *_volStatusMessage[] = {
33  /* 0 */ "ok",
34  /* 1 */ "fault in calling routine",
35  /* 2 */ "out of memory"
36 };
37 /*****************************************************************************/
38 
39 /*****************************************************************************/
45 void volInit(VOL *vol) {
46  if(VOL_TEST) printf("volInit()\n");
47  memset(vol, 0, sizeof(VOL));
49  vol->statmsg=_volStatusMessage[0];
50  vol->orientation=0;
51  vol->dimx=vol->dimy=vol->dimz=0;
52  vol->sizex=vol->sizey=vol->sizez=0;
53  vol->v=(float***)NULL;
54  vol->voxel=(float*)NULL;
55  vol->column=(float*)NULL;
56  vol->row=(float**)NULL;
57  vol->plane=(float***)NULL;
58 }
59 /*****************************************************************************/
60 
61 /*****************************************************************************/
68 void svolInit(SVOL *svol) {
69  if(VOL_TEST) printf("svolInit()\n");
70  memset(svol, 0, sizeof(SVOL));
72  svol->statmsg=_volStatusMessage[0];
73  svol->orientation=0;
74  svol->dimx=svol->dimy=svol->dimz=0;
75  svol->sizex=svol->sizey=svol->sizez=0;
76  svol->scale_factor=1.0;
77  svol->v=(short int***)NULL;
78  svol->voxel=(short int*)NULL;
79  svol->column=(short int*)NULL;
80  svol->row=(short int**)NULL;
81  svol->plane=(short int***)NULL;
82 }
83 /*****************************************************************************/
84 
85 /*****************************************************************************/
91 void volEmpty(VOL *vol) {
92  if(VOL_TEST) printf("volEmpty()\n");
93  if(vol->status<IMG_STATUS_OCCUPIED) return;
94  /* Free up memory */
95  if(vol->_vxl!=NULL) free(vol->_vxl);
96  if(vol->_col!=NULL) free(vol->_col);
97  if(vol->_row!=NULL) free(vol->_row);
98  if(vol->_pln!=NULL) free(vol->_pln);
99  /* Set variables */
100  vol->statmsg=_volStatusMessage[0];
101  vol->orientation=0;
102  vol->dimx=vol->dimy=vol->dimz=0;
103  vol->sizex=vol->sizey=vol->sizez=0;
104  vol->v=(float***)NULL;
105  vol->voxel=(float*)NULL;
106  vol->column=(float*)NULL;
107  vol->row=(float**)NULL;
108  vol->plane=(float***)NULL;
109  /* Set status */
111 }
112 /*****************************************************************************/
113 
114 /*****************************************************************************/
120 void svolEmpty(SVOL *svol) {
121  if(VOL_TEST) printf("svolEmpty()\n");
122  if(svol->status<IMG_STATUS_OCCUPIED) return;
123  /* Free up memory */
124  if(svol->_vxl!=NULL) free(svol->_vxl);
125  if(svol->_col!=NULL) free(svol->_col);
126  if(svol->_row!=NULL) free(svol->_row);
127  if(svol->_pln!=NULL) free(svol->_pln);
128  /* Set variables */
129  svol->statmsg=_volStatusMessage[0];
130  svol->orientation=0;
131  svol->dimx=svol->dimy=svol->dimz=0;
132  svol->sizex=svol->sizey=svol->sizez=0;
133  svol->scale_factor=1.0;
134  svol->v=(short int***)NULL;
135  svol->voxel=(short int*)NULL;
136  svol->column=(short int*)NULL;
137  svol->row=(short int**)NULL;
138  svol->plane=(short int***)NULL;
139  /* Set status */
141 }
142 /*****************************************************************************/
143 
144 /*****************************************************************************/
156 int volAllocate(VOL *vol, int planes, int rows, int columns) {
157  unsigned short int zi, ri;
158  int vxlNr, vi;
159  float **rptr, *cptr;
160 
161  if(VOL_TEST) printf("voiAllocate(*vol, %d, %d, %d)\n", planes, rows, columns);
162  /* Check arguments */
163  if(vol==NULL) return(1); else vol->statmsg=_volStatusMessage[1];
164  if(vol->status==IMG_STATUS_UNINITIALIZED) return(1);
165  if(planes<1 || rows<1 || columns<1) return(2);
166  vxlNr=planes*rows*columns;
167 
168  /* Check if correct volume size is already allocated */
169  if(vol->status>=IMG_STATUS_OCCUPIED) {
170  if(planes==vol->dimz && rows==vol->dimy && columns==vol->dimx) {
171  for(vi=0; vi<vxlNr; vi++) vol->_vxl[vi]=0;
172  return(0); /* that's it */
173  } else {
174  volEmpty(vol);
175  }
176  }
177  /* Allocate memory for volume data */
178  vol->_pln=(float***)malloc(planes*sizeof(float**));
179  if(vol->_pln==NULL) {
180  return(5);}
181  vol->_row=(float**)malloc(planes*rows*sizeof(float*));
182  if(vol->_row==NULL) {
183  free(vol->_pln); return(6);}
184  vol->_col=vol->_vxl=(float*)calloc(planes*rows*columns, sizeof(float));
185  if(vol->_vxl==NULL) {
186  free(vol->_pln); free(vol->_row); return(8);
187  }
188  /* Set data pointers */
189  rptr=vol->_row; cptr=vol->_col;
190  for(zi=0; zi<planes; zi++) {
191  vol->_pln[zi]=rptr;
192  for(ri=0; ri<rows; ri++) {
193  *rptr++=cptr; cptr+=columns;
194  }
195  }
196  vol->v=vol->_pln;
197  vol->plane=vol->_pln;
198  vol->column=vol->_col;
199  vol->row=vol->_row;
200  vol->voxel=vol->_vxl;
201  /* Ok */
202  vol->dimz=planes; vol->dimy=rows; vol->dimx=columns;
203  vol->statmsg=_volStatusMessage[0];
205  return(0);
206 }
207 /*****************************************************************************/
208 
209 /*****************************************************************************/
221 int svolAllocate(SVOL *svol, int planes, int rows, int columns) {
222  unsigned short int zi, ri;
223  int vxlNr, vi;
224  short int **rptr, *cptr;
225 
226  if(VOL_TEST) printf("svoiAllocate(*svol, %d, %d, %d)\n", planes, rows, columns);
227  /* Check arguments */
228  if(svol==NULL) return(1); else svol->statmsg=_volStatusMessage[1];
229  if(svol->status==IMG_STATUS_UNINITIALIZED) return(1);
230  if(planes<1 || rows<1 || columns<1) return(2);
231  vxlNr=planes*rows*columns;
232 
233  /* Check if correct volume size is already allocated */
234  if(svol->status>=IMG_STATUS_OCCUPIED) {
235  if(planes==svol->dimz && rows==svol->dimy && columns==svol->dimx) {
236  for(vi=0; vi<vxlNr; vi++) svol->_vxl[vi]=0;
237  return(0); /* that's it */
238  } else {
239  svolEmpty(svol);
240  }
241  }
242  /* Allocate memory for volume data */
243  svol->_pln=(short int***)malloc(planes*sizeof(short int**));
244  if(svol->_pln==NULL) {
245  return(5);}
246  svol->_row=(short int**)malloc(planes*rows*sizeof(short int*));
247  if(svol->_row==NULL) {
248  free(svol->_pln); return(6);}
249  svol->_col=svol->_vxl=(short int*)calloc(planes*rows*columns, sizeof(short int));
250  if(svol->_vxl==NULL) {
251  free(svol->_pln); free(svol->_row); return(8);
252  }
253  /* Set data pointers */
254  rptr=svol->_row; cptr=svol->_col;
255  for(zi=0; zi<planes; zi++) {
256  svol->_pln[zi]=rptr;
257  for(ri=0; ri<rows; ri++) {
258  *rptr++=cptr; cptr+=columns;
259  }
260  }
261  svol->v=svol->_pln;
262  svol->plane=svol->_pln;
263  svol->column=svol->_col;
264  svol->row=svol->_row;
265  svol->voxel=svol->_vxl;
266  /* Ok */
267  svol->dimz=planes; svol->dimy=rows; svol->dimx=columns;
268  svol->statmsg=_volStatusMessage[0];
270  return(0);
271 }
272 /*****************************************************************************/
273 
274 /*****************************************************************************/
284 int img2vol(IMG *img, VOL *vol, int frame) {
285  int ret;
286  unsigned short int zi, yi, xi, fi;
287 
288  if(VOL_TEST) printf("img2vol(img, %d, vol)\n", frame);
289  /* Check input */
290  if(vol==NULL) return(1);
291  vol->statmsg=_volStatusMessage[1];
292  if(img->status!=IMG_STATUS_OCCUPIED) return(1);
293  if(frame<1 || img->dimt<frame) return(2);
294  if(vol->status==IMG_STATUS_UNINITIALIZED) return(1);
295 
296  /* Allocate memory (if needed) for volume */
297  ret=volAllocate(vol, img->dimz, img->dimy, img->dimx);
298  if(ret) return(ret);
299 
300  /* Copy data */
301  fi=frame-1;
302  vol->orientation=img->orientation;
303  vol->sizex=img->sizex; vol->sizey=img->sizey; vol->sizez=img->sizez;
304  for(zi=0; zi<vol->dimz; zi++)
305  for(yi=0; yi<vol->dimy; yi++)
306  for(xi=0; xi<vol->dimx; xi++)
307  vol->v[zi][yi][xi]=img->m[zi][yi][xi][fi];
308 
309  return(0);
310 }
311 /*****************************************************************************/
312 
313 /*****************************************************************************/
322 int img2svol(IMG *img, SVOL *svol, int frame) {
323  int ret;
324  unsigned short int zi, yi, xi, fi;
325  float fmin, fmax, g;
326 
327  if(VOL_TEST) printf("img2svol(img, %d, svol)\n", frame);
328  /* Check input */
329  if(svol==NULL) return(1);
330  svol->statmsg=_volStatusMessage[1];
331  if(img->status!=IMG_STATUS_OCCUPIED) return(1);
332  if(frame<1 || img->dimt<frame) return(2);
333  if(svol->status==IMG_STATUS_UNINITIALIZED) return(1);
334 
335  /* Allocate memory (if needed) for volume */
336  ret=svolAllocate(svol, img->dimz, img->dimy, img->dimx);
337  if(ret) return(ret);
338 
339  /* Copy data */
340  fi=frame-1;
341  svol->orientation=img->orientation;
342  svol->sizex=img->sizex; svol->sizey=img->sizey; svol->sizez=img->sizez;
343  ret=imgFrameMinMax(img, frame, &fmin, &fmax); if(ret) return(10+ret);
344  if(fabs(fmin)>fabs(fmax)) g=fabs(fmin); else g=fabs(fmax);
345  if(g!=0) g=32766./g; else g=1.0;
346  for(zi=0; zi<svol->dimz; zi++)
347  for(yi=0; yi<svol->dimy; yi++)
348  for(xi=0; xi<svol->dimx; xi++)
349  svol->v[zi][yi][xi]=(short int)temp_roundf(g*img->m[zi][yi][xi][fi]);
350  svol->scale_factor=1.0/g;
351 
352  return(0);
353 }
354 /*****************************************************************************/
355 
356 /*****************************************************************************/
367 int vol2img(VOL *vol, IMG *img, int frame) {
368  unsigned short int zi, yi, xi, fi;
369 
370  if(VOL_TEST) printf("vol2img(vol, img, %d)\n", frame);
371  /* Check input */
372  if(vol==NULL || vol->status!=IMG_STATUS_OCCUPIED) return(1);
373  vol->statmsg=_volStatusMessage[1];
374  if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) return(1);
375  if(frame<1 || img->dimt<frame) return(2);
376  if(img->dimx!=vol->dimx || img->dimy!=vol->dimy) return(3);
377  if(img->dimz!=vol->dimz) return(4);
378 
379  /* Copy data */
380  fi=frame-1;
381  for(zi=0; zi<vol->dimz; zi++)
382  for(yi=0; yi<vol->dimy; yi++)
383  for(xi=0; xi<vol->dimx; xi++)
384  img->m[zi][yi][xi][fi]=vol->v[zi][yi][xi];
385 
386  return(0);
387 }
388 /*****************************************************************************/
389 
390 /*****************************************************************************/
401 int svol2img(SVOL *svol, IMG *img, int frame) {
402  unsigned short int zi, yi, xi, fi;
403 
404  if(VOL_TEST) printf("svol2img(svol, img, %d)\n", frame);
405  /* Check input */
406  if(svol==NULL || svol->status!=IMG_STATUS_OCCUPIED) return(1);
407  svol->statmsg=_volStatusMessage[1];
408  if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) return(1);
409  if(frame<1 || img->dimt<frame) return(2);
410  if(img->dimx!=svol->dimx || img->dimy!=svol->dimy) return(3);
411  if(img->dimz!=svol->dimz) return(4);
412 
413  /* Copy data */
414  fi=frame-1;
415  for(zi=0; zi<svol->dimz; zi++)
416  for(yi=0; yi<svol->dimy; yi++)
417  for(xi=0; xi<svol->dimx; xi++)
418  img->m[zi][yi][xi][fi]=(svol->scale_factor)*(float)svol->v[zi][yi][xi];
419 
420  return(0);
421 }
422 /*****************************************************************************/
423 
424 /*****************************************************************************/
431 void volInfo(VOL *vol, FILE *fp) {
432  if(VOL_TEST) printf("volInfo()\n");
433  if(vol->status<=IMG_STATUS_UNINITIALIZED) {
434  fprintf(fp, "Volume data is not initialized.\n"); return;}
435  if(vol->status==IMG_STATUS_INITIALIZED) {
436  fprintf(fp, "Volume data is initialized but empty.\n"); return;}
437  if(vol->status==IMG_STATUS_ERROR) fprintf(stdout, "Volume data has errors.\n");
438  fprintf(fp, "Volume status: %s\n", vol->statmsg);
439  fprintf(fp, "Patient orientation: %d\n", vol->orientation);
440  fprintf(fp, "Voxel sizes (x, y, z): %g %g %g mm\n",
441  vol->sizex, vol->sizey, vol->sizez);
442  fprintf(fp, "Dimensions (x, y, z): %d %d %d\n",
443  vol->dimx, vol->dimy, vol->dimz);
444  return;
445 }
446 /*****************************************************************************/
447 
448 /*****************************************************************************/
455 void svolInfo(SVOL *svol, FILE *fp) {
456  if(VOL_TEST) printf("svolInfo()\n");
457  if(svol->status<=IMG_STATUS_UNINITIALIZED) {
458  fprintf(fp, "Volume data is not initialized.\n"); return;}
459  if(svol->status==IMG_STATUS_INITIALIZED) {
460  fprintf(fp, "Volume data is initialized but empty.\n"); return;}
461  if(svol->status==IMG_STATUS_ERROR) fprintf(stdout, "Volume data has errors.\n");
462  fprintf(fp, "Volume status: %s\n", svol->statmsg);
463  fprintf(fp, "Patient orientation: %d\n", svol->orientation);
464  fprintf(fp, "Voxel sizes (x, y, z): %g %g %g mm\n",
465  svol->sizex, svol->sizey, svol->sizez);
466  fprintf(fp, "Dimensions (x, y, z): %d %d %d\n",
467  svol->dimx, svol->dimy, svol->dimz);
468  fprintf(fp, "Scale factor: %g\n", svol->scale_factor);
469  return;
470 }
471 /*****************************************************************************/
472 
473 /*****************************************************************************/
480 void volContents(VOL *vol, VOL_RANGE r, FILE *fp) {
481  int zi, yi, xi;
482 
483  if(vol->status!=IMG_STATUS_OCCUPIED) return;
484  if(r.z1<1 || r.y1<1 || r.x1<1) return;
485  if(r.z2<r.z1 || r.y2<r.y1 || r.x2<r.x1) return;
486  if(r.z2>vol->dimz || r.y2>vol->dimy || r.x2>vol->dimx) return;
487 
488  for(zi=r.z1-1; zi<r.z2; zi++) {
489  fprintf(fp, "pl=%03d ", zi+1);
490  for(xi=r.x1-1; xi<r.x2; xi++) fprintf(fp, " x=%05d", xi+1);
491  fprintf(fp, "\n");
492  for(yi=r.y1-1; yi<r.y2; yi++) {
493  fprintf(fp, "y=%05d", yi+1);
494  for(xi=r.x1-1; xi<r.x2; xi++)
495  fprintf(fp, " %7.3f", vol->v[zi][yi][xi]);
496  fprintf(fp, "\n");
497  }
498  }
499 }
500 /*****************************************************************************/
501 
502 /*****************************************************************************/
513 int volMax(VOL *vol, VOL_RANGE r, VOL_PIXEL *p, float *maxv) {
514  int zi, yi, xi;
515 
516  if(vol->status!=IMG_STATUS_OCCUPIED) return(1);
517  if(r.z1<1 || r.y1<1 || r.x1<1) return(2);
518  if(r.z2<r.z1 || r.y2<r.y1 || r.x2<r.x1) return(3);
519  if(r.z2>vol->dimz || r.y2>vol->dimy || r.x2>vol->dimx) return(4);
520 
521  zi=r.z1-1; yi=r.y1-1; xi=r.x1-1; *maxv=vol->v[zi][yi][xi];
522  if(p!=NULL) {p->z=zi+1; p->y=yi+1; p->x=xi+1;}
523  for(zi=r.z1-1; zi<r.z2; zi++) {
524  for(yi=r.y1-1; yi<r.y2; yi++) {
525  for(xi=r.x1-1; xi<r.x2; xi++) if(*maxv<vol->v[zi][yi][xi]) {
526  *maxv=vol->v[zi][yi][xi];
527  if(p!=NULL) {p->z=zi+1; p->y=yi+1; p->x=xi+1;}
528  }
529  }
530  }
531  return(0);
532 }
533 /*****************************************************************************/
534 
535 /*****************************************************************************/
545 int volAvg(VOL *vol, VOL_RANGE r, float *avg) {
546  int zi, yi, xi, n=0;
547 
548  if(vol->status!=IMG_STATUS_OCCUPIED) return(1);
549  if(r.z1<1 || r.y1<1 || r.x1<1) return(2);
550  if(r.z2<r.z1 || r.y2<r.y1 || r.x2<r.x1) return(3);
551  if(r.z2>vol->dimz || r.y2>vol->dimy || r.x2>vol->dimx) return(4);
552 
553  *avg=0.0;
554  for(zi=r.z1-1; zi<r.z2; zi++) {
555  for(yi=r.y1-1; yi<r.y2; yi++) {
556  for(xi=r.x1-1; xi<r.x2; xi++) {
557  *avg+=vol->v[zi][yi][xi]; n++;
558  }
559  }
560  }
561  if(n>0) *avg/=(float)n;
562  return(0);
563 }
564 /*****************************************************************************/
565 
566 /*****************************************************************************/