1 #include "mrilib.h"
2
3 MRI_IMAGE * mri_jointhist( MRI_IMAGE *imp , MRI_IMAGE *imq , byte *mmm ) ;
4
main(int argc,char * argv[])5 int main( int argc , char * argv[] )
6 {
7 int narg , ndset , nvox , nvals,iv ;
8 THD_3dim_dataset *xset , *yset , *mask_dset=NULL , *hset ;
9 byte *mmm=NULL ;
10 MRI_IMAGE *imh , *hdim ;
11 float *har , *hdar ;
12 char *prefix = "jhist" ;
13 THD_ivec3 nxyz ;
14 THD_fvec3 dxyz ;
15
16 /*-- read command line arguments --*/
17
18 if( argc < 3 || strncmp(argv[1],"-help",5) == 0 ){
19 printf("Usage: 3JointHist [options] dset1 dset2\n"
20 "Output = dataset of joint histogram between dset1[0] and\n"
21 " each sub-brick of dset2 (1 histogram per slice).\n"
22 "Options:\n"
23 " -mask mset = Means to use the dataset 'mset' as a mask:\n"
24 " Only voxels with nonzero values in 'mset'\n"
25 " will be averaged from 'dataset'.\n"
26 " -prefix ppp = If you don't know this by now, you shouldn't\n"
27 " even THINK about using this program!\n"
28 ) ;
29 PRINT_COMPILE_DATE ; exit(0) ;
30 }
31
32 narg = 1 ;
33 while( narg < argc && argv[narg][0] == '-' ){
34
35 if( strncmp(argv[narg],"-mask",5) == 0 ){
36 if( mask_dset != NULL ) ERROR_exit("Can't use -mask twice") ;
37 if( narg+1 >= argc ) ERROR_exit("Need argument after -mask") ;
38 mask_dset = THD_open_dataset( argv[++narg] ) ;
39 if( mask_dset == NULL ) ERROR_exit("Can't open -mask %s",argv[narg]) ;
40 narg++ ; continue ;
41 }
42
43 if( strcmp(argv[narg],"-prefix") == 0 ){
44 if( narg+1 >= argc ) ERROR_exit("Need argument after -prefix") ;
45 prefix = argv[++narg] ;
46 narg++ ; continue ;
47 }
48
49 ERROR_exit("Unknown option: %s",argv[narg]) ;
50 }
51
52 /* should have at least 2 more arguments */
53
54 ndset = argc - narg ;
55 if( ndset <= 1 ) ERROR_exit("Need two input datasets") ;
56
57 xset = THD_open_dataset( argv[narg++] ) ;
58 yset = THD_open_dataset( argv[narg++] ) ;
59 if( xset == NULL || yset == NULL )
60 ERROR_exit("Cannot open both input datasets!\n") ;
61
62 if( DSET_NVALS(xset) > 1 )
63 WARNING_message("Will only use sub-brick #0 of 1st input dataset") ;
64 nvals = DSET_NVALS(yset) ;
65
66 nvox = DSET_NVOX(xset) ;
67 if( nvox != DSET_NVOX(yset) )
68 ERROR_exit("Input datasets grid dimensions don't match!\n") ;
69
70 /* make a byte mask from mask dataset */
71
72 if( mask_dset != NULL ){
73 int mcount ;
74 if( DSET_NVOX(mask_dset) != nvox )
75 ERROR_exit("Input and mask datasets are not same dimensions!\n");
76 DSET_load(mask_dset) ; CHECK_LOAD_ERROR(mask_dset) ;
77 mmm = THD_makemask( mask_dset , 0 , 1.0f,-1.0f ) ;
78 mcount = THD_countmask( nvox , mmm ) ;
79 INFO_message("Have %d voxels in the mask\n",mcount) ;
80 if( mcount <= 666 ) ERROR_exit("Mask is too small") ;
81 DSET_delete(mask_dset) ;
82 }
83
84 INFO_message("loading 1st dataset") ;
85 DSET_load(xset) ; CHECK_LOAD_ERROR(xset) ;
86 INFO_message("loading 2nd dataset") ;
87 DSET_load(yset) ; CHECK_LOAD_ERROR(yset) ;
88
89 hset = EDIT_empty_copy(NULL) ;
90 nxyz.ijk[0] = nxyz.ijk[1] = 256 ; nxyz.ijk[2] = nvals ;
91 dxyz.xyz[0] = dxyz.xyz[1] = dxyz.xyz[2] = 1.0f ;
92 EDIT_dset_items( hset ,
93 ADN_prefix , prefix ,
94 ADN_datum_all , MRI_float ,
95 ADN_nvals , 1 ,
96 ADN_type , HEAD_ANAT_TYPE ,
97 ADN_view_type , xset->view_type ,
98 ADN_func_type , ANAT_MRAN_TYPE ,
99 ADN_nxyz , nxyz ,
100 ADN_xyzdel , dxyz ,
101 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
102 ADN_none ) ;
103 EDIT_substitute_brick( hset , 0 , MRI_float , NULL ) ;
104 hdim = DSET_BRICK(hset,0) ; hdar = MRI_FLOAT_PTR(hdim) ;
105
106 fprintf(stderr,"++ histogram-izing") ;
107 for( iv=0 ; iv < nvals ; iv++ ){
108 fprintf(stderr,".") ;
109 imh = mri_jointhist( DSET_BRICK(xset,0), DSET_BRICK(yset,iv), mmm ) ;
110 har = MRI_FLOAT_PTR(imh) ;
111 memcpy( hdar+(256*256*iv) , har , 256*256*sizeof(float) ) ;
112 mri_free(imh) ;
113 }
114 fprintf(stderr,"\n") ;
115
116 DSET_write(hset) ;
117 WROTE_DSET(hset) ;
118 exit(0) ;
119 }
120
121 /*------------------------------------------------------------------------*/
122
mri_jointhist(MRI_IMAGE * imp,MRI_IMAGE * imq,byte * mmm)123 MRI_IMAGE * mri_jointhist( MRI_IMAGE *imp , MRI_IMAGE *imq , byte *mmm )
124 {
125 int nvox , nmmm=0 ;
126 float *rst ;
127 byte *par, *qar ;
128 float fac ;
129 register int ii,jj,kk ;
130 MRI_IMAGE *imqq, *impp , *imh ;
131
132 if( imp == NULL || imq == NULL || imp->nvox != imq->nvox ) return NULL;
133
134 nvox = imp->nvox ;
135
136 impp = (imp->kind==MRI_byte) ? imp : mri_to_byte(imp) ;
137 imqq = (imq->kind==MRI_byte) ? imq : mri_to_byte(imq) ;
138 par = MRI_BYTE_PTR(impp) ;
139 qar = MRI_BYTE_PTR(imqq) ;
140 imh = mri_new( 256,256,MRI_float ) ;
141 rst = MRI_FLOAT_PTR(imh) ;
142
143 if( mmm != NULL ){
144 for( kk=0 ; kk < nvox ; kk++ ) if( mmm[kk] ) nmmm++ ;
145 fac = 1.0f / nmmm ;
146 for( kk=0 ; kk < nvox ; kk++ ){
147 if( mmm[kk] == 0 ) continue ;
148 ii = par[kk] ; jj = qar[kk] ; rst[ii+256*jj] += fac ;
149 }
150 } else {
151 fac = 1.0f / nvox ;
152 for( kk=0 ; kk < nvox ; kk++ ){
153 ii = par[kk] ; jj = qar[kk] ; rst[ii+256*jj] += fac ;
154 }
155 }
156
157 if( impp != imp ) mri_free(impp) ;
158 if( imqq != imq ) mri_free(imqq) ;
159
160 return imh ;
161 }
162