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