1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 /* Cameron Craddock - modification of 3ddot to support linear svm prediction
8 from w (bucket) calculated by 3dsvm */
9
10 #include "mrilib.h"
11 #include <string.h>
12
13 #define SCALE 4000000
14
main(int argc,char * argv[])15 int main( int argc , char * argv[] )
16 {
17 ATR_float * atr_float = NULL;
18 double dxy;
19 int narg = 0;
20 int ndset = 0;
21 int nvox = 0;
22 THD_3dim_dataset *xset = NULL;
23 THD_3dim_dataset *yset = NULL;
24 THD_3dim_dataset *mask_dset= NULL;
25 float mask_bot = 666.0;
26 float mask_top = -666.0 ;
27 double bias_val = 0.0;
28 float *fxar = NULL;
29 float *fyar = NULL;
30 byte *mmm = NULL;
31 int ivx = 0;
32 int ivy = 0;
33 int iv = 0;
34 void *xar = NULL;
35 void *yar = NULL;
36 int itypx = 0;
37 int itypy = 0;
38 int fxar_new = 0;
39 int fyar_new = 0;
40
41 /*-- read command line arguments --*/
42
43 if( argc < 3 || strncmp(argv[1],"-help",5) == 0 )
44 {
45 printf("Usage: 3ddot [options] w dset\n"
46 "Output = linear prediction for w from 3dsvm \n"
47 " - you can use sub-brick selectors on the dsets\n"
48 " - the result is a number printed to stdout\n"
49 "Options:\n"
50 " -mask mset Means to use the dataset 'mset' as a mask:\n"
51 " Only voxels with nonzero values in 'mset'\n"
52 " will be averaged from 'dataset'. Note\n"
53 " that the mask dataset and the input dataset\n"
54 " must have the same number of voxels.\n"
55 );
56
57 printf("\n" MASTER_SHORTHELP_STRING ) ;
58 PRINT_COMPILE_DATE;
59 exit(0);
60 }
61
62 narg = 1;
63 while( narg < argc && argv[narg][0] == '-' )
64 {
65 if( strncmp(argv[narg],"-mask",5) == 0 )
66 {
67 if( mask_dset != NULL )
68 {
69 fprintf(stderr,"*** Cannot have two -mask options!\n") ; exit(1) ;
70 }
71 if( narg+1 >= argc )
72 {
73 fprintf(stderr,"*** -mask option requires a following argument!\n");
74 exit(1) ;
75 }
76 mask_dset = THD_open_dataset( argv[++narg] ) ;
77 if( mask_dset == NULL )
78 {
79 fprintf(stderr,"*** Cannot open mask dataset!\n") ; exit(1) ;
80 }
81 if( DSET_BRICK_TYPE(mask_dset,0) == MRI_complex )
82 {
83 fprintf(stderr,"*** Cannot deal with complex-valued mask dataset!\n");
84 exit(1) ;
85 }
86 narg++;
87 continue;
88 }
89
90 fprintf(stderr,"*** Unknown option: %s\n",argv[narg]);
91 exit(1);
92 }
93
94 ndset = argc - narg ;
95 if( ndset <= 1 )
96 {
97 fprintf(stderr,"*** No input datasets!?\n");
98 exit(1);
99 }
100
101 /* xset is the linear 3dsvm model, yset is the data to be
102 predicted */
103 xset = THD_open_dataset( argv[narg++] );
104 if( xset == NULL )
105 {
106 fprintf(stderr,"*** cannot open first input dataset!\n");
107 DSET_delete(mask_dset) ;
108 exit(1);
109 }
110
111 /* get the B value from the header */
112 atr_float = THD_find_float_atr ( xset->dblk, "3DSVM_B" );
113 if ( atr_float != NULL )
114 {
115 if( atr_float->nfl == 0 )
116 {
117 fprintf(stderr, "Could not find 3DSVM_B. Is the first dataset a bucket output from 3dsvm.\n" );
118 DSET_delete(mask_dset) ;
119 DSET_delete(xset) ;
120 exit(1);
121 }
122 else if( atr_float->nfl > 1 )
123 {
124 fprintf(stderr, "3dsvm_linpredict cannot handle multi-class models.\n" );
125 DSET_delete(mask_dset) ;
126 DSET_delete(xset) ;
127 exit(1);
128 }
129 else
130 {
131 bias_val=atr_float->fl[0];
132 }
133 }
134 else
135 {
136 fprintf(stderr, "Could not find 3DSVM_B. Is the first dataset a bucket output from 3dsvm.\n" );
137 DSET_delete(mask_dset) ;
138 DSET_delete(xset) ;
139 exit(1);
140 }
141
142 yset = THD_open_dataset( argv[narg++] );
143 if( yset == NULL )
144 {
145 fprintf(stderr,"*** cannot open second input datasets!\n");
146 DSET_delete(mask_dset);
147 DSET_delete(xset);
148 exit(1);
149 }
150 if( DSET_NUM_TIMES(xset) > 1 )
151 {
152 fprintf(stderr,"*** cannot use time-dependent datasets!\n");
153 DSET_delete(mask_dset);
154 DSET_delete(xset);
155 DSET_delete(yset);
156 exit(1);
157 }
158 nvox = DSET_NVOX(xset);
159 if( nvox != DSET_NVOX(yset) )
160 {
161 fprintf(stderr,"*** input datasets dimensions don't match!\n");
162 DSET_delete(mask_dset);
163 DSET_delete(xset);
164 DSET_delete(yset);
165 exit(1);
166 }
167 if( !EQUIV_GRIDS(xset,yset) )
168 {
169 WARNING_message("input datasets don't have same grids");
170 }
171
172
173 /* make a byte mask from mask dataset */
174 if( mask_dset != NULL )
175 {
176 int mcount;
177 if( DSET_NVOX(mask_dset) != nvox )
178 {
179 fprintf(stderr,"*** Input and mask datasets are not same dimensions!\n");
180 DSET_delete(mask_dset);
181 DSET_delete(xset);
182 DSET_delete(yset);
183 exit(1);
184 }
185 mmm = THD_makemask( mask_dset, 0, mask_bot, mask_top );
186 mcount = THD_countmask( nvox , mmm );
187 fprintf(stderr,"+++ %d voxels in the mask\n",mcount) ;
188 if( mcount <= 5 )
189 {
190 fprintf(stderr,"*** Mask is too small!\n");
191 DSET_delete(mask_dset);
192 DSET_delete(xset);
193 DSET_delete(yset);
194 exit(1);
195 }
196 DSET_delete(mask_dset) ;
197 }
198
199 /* load bricks */
200
201 DSET_load(xset);
202 CHECK_LOAD_ERROR(xset);
203 ivx = 0 ;
204 itypx = DSET_BRICK_TYPE(xset,ivx) ;
205 xar = DSET_ARRAY(xset,ivx);
206 if( xar == NULL )
207 {
208 fprintf(stderr,"Could not access brick %d in first datset\n", ivx);
209 exit(1);
210 }
211 if( itypx == MRI_float )
212 {
213 fxar = (float *) xar;
214 fxar_new = 0;
215 }
216 else
217 {
218 fxar = (float *) malloc( sizeof(float) * nvox );
219 if( fxar == NULL )
220 {
221 fprintf(stderr,"Could not allocate fxar\n");
222 exit(1);
223 }
224 fxar_new = 1 ;
225 EDIT_coerce_type( nvox, itypx, xar, MRI_float, fxar );
226 }
227
228 DSET_load(yset);
229 CHECK_LOAD_ERROR(yset);
230
231 for( ivy=0; ivy<DSET_NUM_TIMES(yset); ivy++ )
232 {
233 itypy = DSET_BRICK_TYPE(yset,ivy);
234 yar = DSET_ARRAY(yset,ivy);
235 if( yar == NULL )
236 {
237 fprintf(stderr,"Could not access brick %d in second datset\n", ivy);
238 if( fxar_new == 1 )
239 free(fxar);
240 fxar=NULL;
241 DSET_delete(mask_dset);
242 DSET_delete(xset);
243 DSET_delete(yset);
244 exit(1);
245 }
246
247 if( itypy == MRI_float )
248 {
249 fyar = (float *) yar ; fyar_new = 0 ;
250 }
251 else
252 {
253 if( fyar == NULL )
254 {
255 fyar = (float *) malloc( sizeof(float) * nvox );
256 if( fyar == NULL )
257 {
258 fprintf(stderr,"Could not allocate fyar\n");
259 if( fxar_new == 1 )
260 free( fxar);
261 fyar=NULL;
262 DSET_delete(mask_dset);
263 DSET_delete(xset);
264 DSET_delete(yset);
265 exit(1);
266 }
267 fyar_new = 1;
268 }
269 EDIT_coerce_type( nvox, itypy, yar, MRI_float, fyar ) ;
270 }
271
272 dxy = 0.0;
273 for( iv = 0; iv<nvox; iv++ )
274 {
275 if(( mmm == NULL ) || ( mmm[iv] == 1 ))
276 {
277 dxy += ((double)fxar[iv]/(double)SCALE)*((double)fyar[iv]);
278 }
279 }
280 dxy -= bias_val;
281 printf( "%g\n",dxy );
282 }
283
284 if ( fxar_new == 1 ) free( fxar );
285 if ( fyar_new == 1 ) free( fyar );
286 if ( mmm != NULL ) free ( mmm );
287 DSET_delete(mask_dset);
288 DSET_delete(xset);
289 DSET_delete(yset);
290 exit(0) ;
291 }
292