1 /*
2 Description
3 */
4
5
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <math.h>
9 #include <unistd.h>
10 #include <debugtrace.h>
11 #include "mrilib.h"
12 #include "3ddata.h"
13 //#include "rsfc.h"
14 //#include <gsl/gsl_rng.h>
15 #include "DoTrackit.h"
16
usage_FUNCNAME(int detail)17 void usage_FUNCNAME(int detail)
18 {
19 printf(
20 "\n"
21 " *** \n"
22 " *** \n"
23 "\n"
24 "* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\n"
25 " \n"
26 " + USAGE: ***\n"
27 "\n"
28 "* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\n"
29 "\n"
30 " + COMMAND: *** \n"
31 "\n"
32 "* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\n"
33 "\n"
34 " + RUNNING, need to provide:\n"
35 " *** \n"
36 "\n"
37 "* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\n"
38 " + OUTPUT: \n"
39 " *** \n"
40 "\n"
41 "* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\n"
42 "\n"
43 " + EXAMPLE:\n"
44 " *** \n"
45 "\n"
46 "* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\n"
47 "\n"
48 " reference *** \n"
49 "____________________________________________________________________________\n"
50 );
51 return;
52 }
53
main(int argc,char * argv[])54 int main(int argc, char *argv[]) {
55 int i,j,k,m,n,mm;
56 int idx=0;
57 int iarg;
58
59 THD_3dim_dataset *insetA = NULL;
60 THD_3dim_dataset *MASK=NULL;
61 char *prefix="PREFIX" ;
62 // char in_name[THD_MAX_NAME];
63
64 // FILE *fout0, *fout1;
65
66 int Nvox=-1; // tot number vox
67 int *Dim=NULL;
68 int ***mskd=NULL; // define mask of where time series are nonzero
69
70 int TEST_OK = 0;
71
72 mainENTRY("3dFUNCNAME"); machdep();
73
74 // ****************************************************************
75 // ****************************************************************
76 // load AFNI stuff
77 // ****************************************************************
78 // ****************************************************************
79
80 // INFO_message("version: NU");
81
82 /** scan args **/
83 if (argc == 1) { usage_FUNCNAME(1); exit(0); }
84 iarg = 1;
85 while( iarg < argc && argv[iarg][0] == '-' ){
86 if( strcmp(argv[iarg],"-help") == 0 ||
87 strcmp(argv[iarg],"-h") == 0 ) {
88 usage_FUNCNAME(strlen(argv[iarg])>3 ? 2:1);
89 exit(0);
90 }
91
92 // NO ARG:
93 if( strcmp(argv[iarg],"-TESTING") == 0) {
94 TEST_OK=1;
95 iarg++ ; continue ;
96 }
97
98 if( strcmp(argv[iarg],"-prefix") == 0 ){
99 iarg++ ; if( iarg >= argc )
100 ERROR_exit("Need argument after '-prefix'");
101 prefix = strdup(argv[iarg]) ;
102 if( !THD_filename_ok(prefix) )
103 ERROR_exit("Illegal name after '-prefix'");
104 iarg++ ; continue ;
105 }
106
107 if( strcmp(argv[iarg],"-insetA") == 0 ){
108 iarg++ ; if( iarg >= argc )
109 ERROR_exit("Need argument after '-insetA'");
110
111 insetA = THD_open_dataset(argv[iarg]);
112 if( (insetA == NULL ))
113 ERROR_exit("Can't open time series dataset '%s'.",
114 argv[iarg]);
115
116 Dim = (int *)calloc(4,sizeof(int));
117 DSET_load(insetA); CHECK_LOAD_ERROR(insetA);
118 Nvox = DSET_NVOX(insetA) ;
119 Dim[0] = DSET_NX(insetA); Dim[1] = DSET_NY(insetA);
120 Dim[2] = DSET_NZ(insetA); Dim[3] = DSET_NVALS(insetA);
121
122 iarg++ ; continue ;
123 }
124
125
126 if( strcmp(argv[iarg],"-mask") == 0 ){
127 iarg++ ; if( iarg >= argc )
128 ERROR_exit("Need argument after '-mask'");
129
130 MASK = THD_open_dataset(argv[iarg]);
131 if( MASK == NULL )
132 ERROR_exit("Can't open time series dataset '%s'.",
133 argv[iarg]);
134
135 DSET_load(MASK); CHECK_LOAD_ERROR(MASK);
136
137 iarg++ ; continue ;
138 }
139
140 // EXAMPLE: option with numerical input: ATOF
141 /*if( strcmp(argv[iarg],"-neigh_Y") == 0 ){
142 iarg++ ; if( iarg >= argc )
143 ERROR_exit("Need argument after '-nneigh'");
144
145 NEIGH_R[1] = atof(argv[iarg]);
146
147 iarg++ ; continue ;
148 }*/
149
150
151 ERROR_message("Bad option '%s'\n",argv[iarg]) ;
152 suggest_best_prog_option(argv[0], argv[iarg]);
153 exit(1);
154 }
155
156
157 // TEST BASIC INPUT PROPERTIES
158 if (iarg < 3) {
159 ERROR_message("Too few options. Try -help for details.\n");
160 exit(1);
161 }
162
163 if( !TEST_OK ) {
164 ERROR_message("HEY! Just testing/building mode right now!\n");
165 exit(5);
166 }
167
168
169
170
171 // ****************************************************************
172 // ****************************************************************
173 // pre-stuff, make storage
174 // ****************************************************************
175 // ****************************************************************
176
177 mskd = (int ***) calloc( Dim[0], sizeof(int **) );
178 for ( i = 0 ; i < Dim[0] ; i++ )
179 mskd[i] = (int **) calloc( Dim[1], sizeof(int *) );
180 for ( i = 0 ; i < Dim[0] ; i++ )
181 for ( j = 0 ; j < Dim[1] ; j++ )
182 mskd[i][j] = (int *) calloc( Dim[2], sizeof(int) );
183
184
185
186 // *************************************************************
187 // *************************************************************
188 // Beginning of main loops
189 // *************************************************************
190 // *************************************************************
191
192
193 // go through once: define data vox
194 idx = 0;
195 for( k=0 ; k<Dim[2] ; k++ )
196 for( j=0 ; j<Dim[1] ; j++ )
197 for( i=0 ; i<Dim[0] ; i++ ) {
198 if( MASK ) {
199 if( THD_get_voxel(MASK,idx,0)>0 )
200 mskd[i][j][k] = 1;
201 }
202 else
203 if( fabs(THD_get_voxel(insetA,idx,0))+
204 fabs(THD_get_voxel(insetA,idx,1))+
205 fabs(THD_get_voxel(insetA,idx,2))+
206 fabs(THD_get_voxel(insetA,idx,3))+
207 fabs(THD_get_voxel(insetA,idx,4)) > EPS_V)
208 mskd[i][j][k] = 1;
209 idx+= 1; // skip, and mskd and KW are both still 0 from calloc
210 }
211
212
213 // **************************************************************
214 // **************************************************************
215 // Store and output
216 // **************************************************************
217 // **************************************************************
218
219
220 // ************************************************************
221 // ************************************************************
222 // Freeing
223 // ************************************************************
224 // ************************************************************
225
226 if(insetA){
227 DSET_delete(insetA);
228 free(insetA);
229 }
230
231 if( MASK ) {
232 DSET_delete(MASK);
233 free(MASK);
234 }
235
236
237 if(mskd) {
238 for( i=0 ; i<Dim[0] ; i++)
239 for( j=0 ; j<Dim[1] ; j++) {
240 free(mskd[i][j]);
241 }
242 for( i=0 ; i<Dim[0] ; i++) {
243 free(mskd[i]);
244 }
245 free(mskd);
246 }
247
248 if(prefix)
249 free(prefix);
250
251 if(Dim)
252 free(Dim);
253
254 return 0;
255 }
256