1 #include "mrilib.h"
2
3 static char prefix[THD_MAX_PREFIX] = "tnorm" ;
4 static int datum = MRI_float ;
5 static int mode = 2 ;
6 static int polort = -1 ;
7 static int dmode = 2 ;
8
9 static void NORM_tsfunc( double tzero , double tdelta ,
10 int npts , float ts[] , double ts_mean ,
11 double ts_slope , void *ud , int nbriks, float *val ) ;
12
main(int argc,char * argv[])13 int main( int argc , char *argv[] )
14 {
15 THD_3dim_dataset *old_dset , *new_dset ; /* input and output datasets */
16 int nopt, ii , nvals ;
17
18 /*----- Help the blessedly ignorant user? -----*/
19
20 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
21 printf("Usage: 3dTnorm [options] dataset\n"
22 "Takes each voxel time series and normalizes it\n"
23 "(by multiplicative scaling) -- in some sense.\n"
24 "\n"
25 "Options:\n"
26 " -prefix p = use string 'p' for the prefix of the\n"
27 " output dataset [DEFAULT = 'tnorm']\n"
28 " -norm2 = L2 normalize (sum of squares = 1) [DEFAULT]\n"
29 " -normR = normalize so sum of squares = number of time points\n"
30 " * e.g., so RMS = 1.\n"
31 " -norm1 = L1 normalize (sum of absolute values = 1)\n"
32 " -normx = Scale so max absolute value = 1 (L_infinity norm)\n"
33 " -polort p = Detrend with polynomials of order p before normalizing\n"
34 " [DEFAULT = don't do this]\n"
35 " * Use '-polort 0' to remove the mean, for example\n"
36 " -L1fit = Detrend with L1 regression (L2 is the default)\n"
37 " * This option is here just for the hell of it\n"
38 "\n"
39 "Notes:\n"
40 "* Each voxel is processed separately\n"
41 "* A voxel that is all zero will be unchanged (duh)\n"
42 "* Output dataset is in float format, no matter what the input format\n"
43 "* This program is for producing regressors to use in 3dTfitter\n"
44 "* Also see programs 1dnorm and 3dcalc\n"
45 ) ;
46 PRINT_COMPILE_DATE ; exit(0) ;
47 }
48
49 /* bureaucracy */
50
51 mainENTRY("3dTnorm main"); machdep(); AFNI_logger("3dTnorm",argc,argv);
52 PRINT_VERSION("3dTnorm"); AUTHOR("RW Cox");
53
54 /*--- scan command line for options ---*/
55
56 nopt = 1 ;
57 while( nopt < argc && argv[nopt][0] == '-' ){
58
59 /*-- prefix --*/
60
61 if( strcmp(argv[nopt],"-prefix") == 0 ){
62 if( ++nopt >= argc ) ERROR_exit("-prefix needs an argument!");
63 MCW_strncpy(prefix,argv[nopt],THD_MAX_PREFIX) ;
64 if( !THD_filename_ok(prefix) )
65 ERROR_exit("%s is not a valid prefix!",prefix);
66 nopt++ ; continue ;
67 }
68
69 if( strncmp(argv[nopt],"-norm",5) == 0 ){
70 if( argv[nopt][5] == '1' ) mode = 1 ;
71 else if( argv[nopt][5] == 'x' ) mode = 666 ;
72 else if( argv[nopt][5] == '2' ) mode = 2 ;
73 else if( argv[nopt][5] == 'R' ) mode = 3 ;
74 else ERROR_message("Don't understand option '%s'",argv[nopt]) ;
75 nopt++ ; continue ;
76 }
77
78 if( strcmp(argv[nopt],"-polort") == 0 ){
79 if( ++nopt >= argc ) ERROR_exit("-polort needs an argument!");
80 polort = (int)strtod(argv[nopt],NULL) ;
81 nopt++ ; continue ;
82 }
83
84 if( strcasecmp(argv[nopt],"-L1fit") == 0 ){
85 dmode = 1 ; nopt++ ; continue ;
86 }
87 if( strcasecmp(argv[nopt],"-L2fit") == 0 ){
88 dmode = 2 ; nopt++ ; continue ;
89 }
90
91 /*-- Quien sabe'? --*/
92
93 ERROR_message("Unknown option: %s",argv[nopt]) ;
94 suggest_best_prog_option(argv[0],argv[nopt]);
95 exit(1);
96 }
97
98 /*----- read input dataset -----*/
99
100 if( nopt >= argc ) ERROR_exit(" No input dataset!?") ;
101
102 old_dset = THD_open_dataset( argv[nopt] ) ;
103 if( !ISVALID_DSET(old_dset) )
104 ERROR_exit("Can't open dataset %s",argv[nopt]);
105 DSET_load(old_dset) ;
106 if( !DSET_LOADED(old_dset) )
107 ERROR_exit("Can't load dataset %s",argv[nopt]) ;
108
109 nvals = DSET_NVALS(old_dset) ;
110 if( nvals < 2 || nvals <= polort )
111 ERROR_exit("Need at least %d time points, but have only %d",
112 MAX(2,polort+1) , nvals ) ;
113
114 nopt++ ;
115 if( nopt < argc )
116 WARNING_message("Trailing inputs on command line ignored: %s ...",argv[nopt]) ;
117
118 /*------------- ready to compute new dataset -----------*/
119
120 new_dset = MAKER_4D_to_typed_fbuc(
121 old_dset , /* input dataset */
122 prefix , /* output prefix */
123 datum , /* output datum */
124 0 , /* ignore count */
125 0 , /* don't detrend */
126 nvals , /* number of briks in output */
127 NORM_tsfunc , /* timeseries processor */
128 NULL, /* data for tsfunc */
129 NULL, /* mask */
130 0 /* Allow auto scaling of output */
131 ) ;
132
133 if( new_dset != NULL ){
134 tross_Copy_History( old_dset , new_dset ) ;
135 tross_Make_History( "3dTnorm" , argc,argv , new_dset ) ;
136 EDIT_dset_items( new_dset ,
137 ADN_ntt , nvals ,
138 ADN_ttorg , DSET_TIMEORIGIN(old_dset) ,
139 ADN_ttdel , DSET_TR(old_dset) ,
140 ADN_tunits , UNITS_SEC_TYPE ,
141 NULL ) ;
142 DSET_write(new_dset) ; WROTE_DSET(new_dset) ;
143 } else {
144 ERROR_exit("Unable to compute output dataset for some unknowable reason :-(") ;
145 }
146
147 exit(0) ;
148 }
149
150 /**********************************************************************
151 Function that does the real work
152 ***********************************************************************/
153
NORM_tsfunc(double tzero,double tdelta,int npts,float ts[],double ts_mean,double ts_slope,void * ud,int nbriks,float * val)154 static void NORM_tsfunc( double tzero, double tdelta ,
155 int npts, float ts[],
156 double ts_mean, double ts_slope,
157 void *ud, int nbriks, float *val )
158 {
159 int ii , nval ;
160
161 /** is this a "notification"? **/
162
163 if( val == NULL ){
164 #if 0
165 if( npts > 0 ){ /* the "start notification" */
166 } else { /* the "end notification" */
167 }
168 #endif
169 return ;
170 }
171
172 /** do the work **/
173
174 nval = MIN(nbriks,npts) ;
175 memcpy( val , ts , sizeof(float)*nval ) ;
176
177 if( polort >= 0 ){
178 switch( dmode ){
179 default:
180 case 2: THD_generic_detrend_LSQ( nval,val, polort, 0,NULL,NULL ); break;
181 case 1: THD_generic_detrend_L1 ( nval,val, polort, 0,NULL,NULL ); break;
182 }
183 }
184
185 switch( mode ){
186 default:
187 case 2: THD_normalize( nval,val ) ; break ;
188 case 1: THD_normL1 ( nval,val ) ; break ;
189 case 3: THD_normRMS ( nval,val ) ; break ;
190 case 666: THD_normmax ( nval,val ) ; break ;
191 }
192
193 return ;
194 }
195