My Project
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/*****************************************************************************/
33 /* 0 */ "ok",
34 /* 1 */ "fault in calling routine",
35 /* 2 */ "out of memory"
36};
37/*****************************************************************************/
38
39/*****************************************************************************/
45void volInit(VOL *vol) {
46 if(VOL_TEST) printf("volInit()\n");
47 memset(vol, 0, sizeof(VOL));
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/*****************************************************************************/
68void svolInit(SVOL *svol) {
69 if(VOL_TEST) printf("svolInit()\n");
70 memset(svol, 0, sizeof(SVOL));
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/*****************************************************************************/
91void 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 */
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/*****************************************************************************/
120void 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/*****************************************************************************/
156int 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;
205 return(0);
206}
207/*****************************************************************************/
208
209/*****************************************************************************/
221int 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/*****************************************************************************/
284int 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);
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/*****************************************************************************/
322int 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/*****************************************************************************/
367int 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);
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/*****************************************************************************/
401int 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/*****************************************************************************/
431void volInfo(VOL *vol, FILE *fp) {
432 if(VOL_TEST) printf("volInfo()\n");
434 fprintf(fp, "Volume data is not initialized.\n"); return;}
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/*****************************************************************************/
455void svolInfo(SVOL *svol, FILE *fp) {
456 if(VOL_TEST) printf("svolInfo()\n");
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/*****************************************************************************/
480void 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/*****************************************************************************/
513int 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/*****************************************************************************/
545int 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/*****************************************************************************/
#define IMG_STATUS_OCCUPIED
Definition img.h:73
#define IMG_STATUS_INITIALIZED
Definition img.h:72
#define IMG_STATUS_UNINITIALIZED
Definition img.h:71
#define IMG_STATUS_ERROR
Definition img.h:74
int imgFrameMinMax(IMG *img, int frame, float *minvalue, float *maxvalue)
Definition imgmax.c:147
Definition img.h:156
float sizex
Definition img.h:208
unsigned short int dimx
Definition img.h:261
float **** m
Definition img.h:274
char status
Definition img.h:164
float sizey
Definition img.h:210
unsigned short int dimz
Definition img.h:265
unsigned short int dimy
Definition img.h:263
int orientation
Definition img.h:188
float sizez
Definition img.h:212
Definition vol.h:62
float sizey
Definition vol.h:71
unsigned short int dimx
Definition vol.h:73
short int * voxel
Definition vol.h:79
short int * _col
Definition vol.h:77
short int ** row
Definition vol.h:79
unsigned short int dimz
Definition vol.h:73
int orientation
Definition vol.h:69
short int * _vxl
Definition vol.h:77
short int *** _pln
Definition vol.h:77
float sizez
Definition vol.h:71
unsigned short int dimy
Definition vol.h:73
short int ** _row
Definition vol.h:77
short int *** plane
Definition vol.h:79
float scale_factor
Definition vol.h:75
char status
Definition vol.h:65
short int * column
Definition vol.h:79
char * statmsg
Definition vol.h:67
short int *** v
Definition vol.h:79
float sizex
Definition vol.h:71
int x
Definition vol.h:25
int z
Definition vol.h:27
int y
Definition vol.h:26
int x2
Definition vol.h:30
int z1
Definition vol.h:32
int y1
Definition vol.h:31
int y2
Definition vol.h:31
int z2
Definition vol.h:32
int x1
Definition vol.h:30
Definition vol.h:40
char status
Definition vol.h:43
float * _vxl
Definition vol.h:53
float * column
Definition vol.h:55
unsigned short int dimx
Definition vol.h:51
unsigned short int dimy
Definition vol.h:51
float *** plane
Definition vol.h:55
int orientation
Definition vol.h:47
float *** _pln
Definition vol.h:53
float ** _row
Definition vol.h:53
unsigned short int dimz
Definition vol.h:51
float * voxel
Definition vol.h:55
float sizex
Definition vol.h:49
char * statmsg
Definition vol.h:45
float ** row
Definition vol.h:55
float *** v
Definition vol.h:55
float sizez
Definition vol.h:49
float sizey
Definition vol.h:49
float * _col
Definition vol.h:53
char * _volStatusMessage[]
Definition vol.c:32
int volMax(VOL *vol, VOL_RANGE r, VOL_PIXEL *p, float *maxv)
Definition vol.c:513
void svolInit(SVOL *svol)
Definition vol.c:68
void volContents(VOL *vol, VOL_RANGE r, FILE *fp)
Definition vol.c:480
int svolAllocate(SVOL *svol, int planes, int rows, int columns)
Definition vol.c:221
void volInit(VOL *vol)
Definition vol.c:45
void svolEmpty(SVOL *svol)
Definition vol.c:120
int svol2img(SVOL *svol, IMG *img, int frame)
Definition vol.c:401
int volAllocate(VOL *vol, int planes, int rows, int columns)
Definition vol.c:156
void volInfo(VOL *vol, FILE *fp)
Definition vol.c:431
void volEmpty(VOL *vol)
Definition vol.c:91
void svolInfo(SVOL *svol, FILE *fp)
Definition vol.c:455
int img2vol(IMG *img, VOL *vol, int frame)
Definition vol.c:284
int img2svol(IMG *img, SVOL *svol, int frame)
Definition vol.c:322
int volAvg(VOL *vol, VOL_RANGE r, float *avg)
Definition vol.c:545
int vol2img(VOL *vol, IMG *img, int frame)
Definition vol.c:367
int VOL_TEST
Definition vol.h:20