1 #include "mrilib.h"
2
THD_mean_dataset(int nds,THD_3dim_dataset ** dsin,int ivbot,int ivtop,int verb)3 THD_3dim_dataset * THD_mean_dataset( int nds, THD_3dim_dataset **dsin, int ivbot, int ivtop, int verb )
4 {
5 THD_3dim_dataset *dsout ;
6 int nbr,dd , nx,ny,nz,nxyz , vv,kk , ndd , ii ;
7 float **br , *bb , fac ;
8
9 ENTRY("THD_mean_dataset") ;
10
11 /* check inputs */
12
13 if( nds <= 0 || dsin == NULL || !ISVALID_DSET(dsin[0]) ) RETURN(NULL) ;
14 if( ivbot < 0 ) ivbot = 0 ;
15 if( ivtop < ivbot || ivtop >= DSET_NVALS(dsin[0]) ) ivtop = DSET_NVALS(dsin[0])-1 ;
16 if( ivtop < ivbot ) RETURN(NULL) ;
17
18 nbr = ivtop-ivbot+1 ;
19
20 nx = DSET_NX(dsin[0]) ;
21 ny = DSET_NY(dsin[0]) ;
22 nz = DSET_NZ(dsin[0]) ; nxyz = nx*ny*nz ;
23
24 for( dd=0 ; dd < nds ; dd++ ){
25 if( !ISVALID_DSET(dsin[dd]) ) RETURN(NULL) ;
26 if( ivtop >= DSET_NVALS(dsin[dd]) ) RETURN(NULL) ;
27 if( DSET_NX(dsin[dd]) != nx ||
28 DSET_NY(dsin[dd]) != ny || DSET_NZ(dsin[dd]) != nz ) RETURN(NULL) ;
29 }
30
31 /* create output dataset */
32
33 dsout = EDIT_empty_copy(dsin[0]) ;
34 EDIT_dset_items( dsout ,
35 ADN_nvals , nbr ,
36 ADN_brick_fac , NULL ,
37 ADN_none ) ;
38 br = (float **)malloc(sizeof(float *)*nbr) ;
39 for( vv=0 ; vv < nbr ; vv++ ){
40 br[vv] = (float *)calloc(sizeof(float),nxyz) ;
41 EDIT_substitute_brick( dsout , vv , MRI_float , br[vv] ) ;
42 }
43
44 /* loop over input datasets */
45
46 for( ndd=dd=0 ; dd < nds ; dd++ ){
47 DSET_load(dsin[dd]) ;
48 if( !DSET_LOADED(dsin[dd]) ){
49 ERROR_message("Can't load dataset %s in THD_mean_dataset()",DSET_BRIKNAME(dsin[dd])) ;
50 continue ;
51 }
52 if( verb ) fprintf(stderr,".") ;
53
54 for( vv=0 ; vv < nbr ; vv++ ){
55 kk = vv+ivbot ; bb = br[vv] ;
56 switch( DSET_BRICK_TYPE(dsin[dd],kk) ){
57 default:
58 ERROR_message("THD_mean_dataset(): can't use sub-brick with datum=%d", DSET_BRICK_TYPE(dsin[dd],kk));
59 ndd-- ;
60 break ;
61
62 case MRI_float:{
63 float *pp = (float *) DSET_ARRAY(dsin[dd],kk) ;
64 fac = DSET_BRICK_FACTOR(dsin[dd],kk) ; if( fac == 0.0f ) fac = 1.0f ;
65 for( ii=0 ; ii < nxyz ; ii++ ) bb[ii] += fac * pp[ii] ;
66 }
67 break ;
68
69 case MRI_short:{
70 short *pp = (short *) DSET_ARRAY(dsin[dd],kk) ;
71 fac = DSET_BRICK_FACTOR(dsin[dd],kk) ; if( fac == 0.0f ) fac = 1.0f ;
72 for( ii=0 ; ii < nxyz ; ii++ ) bb[ii] += fac * pp[ii] ;
73 }
74 break ;
75
76 case MRI_byte:{
77 byte *pp = (byte *) DSET_ARRAY(dsin[dd],kk) ;
78 fac = DSET_BRICK_FACTOR(dsin[dd],kk) ; if( fac == 0.0f ) fac = 1.0f ;
79 for( ii=0 ; ii < nxyz ; ii++ ) bb[ii] += fac * pp[ii] ;
80 }
81 break ;
82 } /* end of switch on sub-brick type */
83 } /* end of sub-brick loop */
84
85 ndd++ ; DSET_unload(dsin[dd]) ; /* one more done */
86 }
87
88 if( ndd == 0 ) RETURN(dsout) ;
89
90 /* scale to be mean */
91
92 fac = 1.0f / ndd ;
93 for( vv=0 ; vv < nbr ; vv++ ){
94 bb = br[vv] ;
95 for( ii=0 ; ii < nxyz ; ii++ ) bb[ii] *= fac ;
96 }
97
98 free(br) ; RETURN(dsout) ;
99 }
100