1 /*
2 Program to calculate ALFF, mALFF, fALFF (as well as RSFA and similar
3 variations) for resting state time series. This program is
4 **heavily** based on the existing 3dBandPass by RW Cox, with the
5 amendments to calculate RSFC parameters written by PA Taylor (July,
6 2012).
7
8 All options of 3dBandPass may be used here (with a couple other
9 parameter options, as well): essentially, the motivation of this
10 program is to produce ALFF, etc. values of the actual RSFC time
11 series that you calculate. Therefore, all the 3dBandPass processing
12 you normally do en route to making your final `resting state time
13 series' is done here to generate your LFFs, from which the
14 amplitudes in the LFF band are calculated at the end. In order to
15 calculate fALFF, the same initial time series are put through the
16 same processing steps which you have chosen but *without* the
17 bandpass part; the spectrum of this second time series is used to
18 calculate the fALFF denominator.
19
20 Sept. 2012: improving some memory stuff
21 Aug. 2016: change >f_N initialization to avoid division errors when
22 there are *a lot* of time points.
23
24 */
25
26
27
28 #include "mrilib.h"
29 #include <gsl/gsl_errno.h> // @@ for ffting here
30 #include <gsl/gsl_fft_real.h>
31 #include <gsl/gsl_fft_halfcomplex.h>
32
33
34 #ifdef USE_OMP
35 #include "mri_blur3d_variable.c"
36 #include "thd_satcheck.c"
37 #endif
38
39 void THD_vectim_localpv( MRI_vectim *mrv , float rad ) ;
40 int THD_vectim_subset_pv( MRI_vectim *mrv, int nind, int *ind,
41 float *ar, unsigned short xran[] ) ;
42
43 /*----------------------------------------------------------------------------*/
44
main(int argc,char * argv[])45 int main( int argc , char * argv[] )
46 {
47 int do_norm=0 , qdet=2 , have_freq=0 , do_automask=0 ;
48 float dt=0.0f , fbot=0.0f,ftop=999999.9f , blur=0.0f ;
49 MRI_IMARR *ortar=NULL ; MRI_IMAGE *ortim=NULL ;
50 THD_3dim_dataset **ortset=NULL ; int nortset=0 ;
51 THD_3dim_dataset *inset=NULL , *outset=NULL;
52 char *prefix="RSFC" ;
53 byte *mask=NULL ;
54 int mask_nx=0,mask_ny=0,mask_nz=0,nmask , verb=1 ,
55 nx,ny,nz,nvox , nfft=0 , kk ;
56 float **vec , **ort=NULL ; int nort=0 , vv , nopt , ntime ;
57 MRI_vectim *mrv ;
58 float pvrad=0.0f ; int nosat=0 ;
59 int do_despike=0 ;
60
61 // @@ non-BP variables
62 float fbotALL=0.0f, ftopALL=999999.9f; // do full range version
63 int NumDen = 0; // switch for doing numerator or denom
64 THD_3dim_dataset *outsetALL=NULL ;
65 int m, mm;
66 float delf; // harmonics
67 int ind_low,ind_high,N_ny, ctr;
68 float sqnt,nt_fac;
69 gsl_fft_real_wavetable *real1, *real2; // GSL stuff
70 gsl_fft_real_workspace *work;
71 double *series1, *series2;
72 double *xx1,*xx2;
73 float numer,denom,val;
74 float *alff=NULL,*malff=NULL,*falff=NULL,
75 *rsfa=NULL,*mrsfa=NULL,*frsfa=NULL; // values
76 float meanALFF=0.0f,meanRSFA=0.0f; // will be for mean in brain region
77 THD_3dim_dataset *outsetALFF=NULL;
78 THD_3dim_dataset *outsetmALFF=NULL;
79 THD_3dim_dataset *outsetfALFF=NULL;
80 THD_3dim_dataset *outsetRSFA=NULL;
81 THD_3dim_dataset *outsetmRSFA=NULL;
82 THD_3dim_dataset *outsetfRSFA=NULL;
83 char out_lff[THD_MAX_NAME];
84 char out_alff[THD_MAX_NAME];
85 char out_malff[THD_MAX_NAME];
86 char out_falff[THD_MAX_NAME];
87 char out_rsfa[THD_MAX_NAME];
88 char out_mrsfa[THD_MAX_NAME];
89 char out_frsfa[THD_MAX_NAME];
90 char out_unBP[THD_MAX_NAME];
91 int SERIES_OUT = 1;
92 int UNBP_OUT = 0;
93 int DO_RSFA = 1;
94 int BP_LAST = 0; // option for only doing filter to LFFs at very end of proc
95 float de_rsfa=0.0f,nu_rsfa=0.0f;
96 double pow1=0.0,pow2=0.0;
97
98 /*-- help? --*/
99
100 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
101 printf(
102 "\n Program to calculate common resting state functional connectivity (RSFC)\n"
103 " parameters (ALFF, mALFF, fALFF, RSFA, etc.) for resting state time\n"
104 " series. This program is **heavily** based on the existing\n"
105 " 3dBandPass by RW Cox, with the amendments to calculate RSFC\n"
106 " parameters written by PA Taylor (July, 2012).\n"
107 " This program is part of FATCAT (Taylor & Saad, 2013) in AFNI. Importantly,\n"
108 " its functionality can be included in the `afni_proc.py' processing-script \n"
109 " generator; see that program's help file for an example including RSFC\n"
110 " and spectral parameter calculation via the `-regress_RSFC' option.\n"
111 "\n"
112 "* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\n"
113 "\n"
114 " All options of 3dBandPass may be used here (with a couple other\n"
115 " parameter options, as well): essentially, the motivation of this\n"
116 " program is to produce ALFF, etc. values of the actual RSFC time\n"
117 " series that you calculate. Therefore, all the 3dBandPass processing\n"
118 " you normally do en route to making your final `resting state time\n"
119 " series' is done here to generate your LFFs, from which the\n"
120 " amplitudes in the LFF band are calculated at the end. In order to\n"
121 " calculate fALFF, the same initial time series are put through the\n"
122 " same processing steps which you have chosen but *without* the\n"
123 " bandpass part; the spectrum of this second time series is used to\n"
124 " calculate the fALFF denominator.\n"
125 " \n"
126 " For more information about each RSFC parameter, see, e.g.: \n"
127 " ALFF/mALFF -- Zang et al. (2007),\n"
128 " fALFF -- Zou et al. (2008),\n"
129 " RSFA -- Kannurpatti & Biswal (2008).\n"
130 "\n"
131 "* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\n"
132 "\n"
133 " + USAGE: 3dRSFC [options] fbot ftop dataset\n"
134 "\n"
135 "* One function of this program is to prepare datasets for input\n"
136 " to 3dSetupGroupInCorr. Other uses are left to your imagination.\n"
137 "\n"
138 "* 'dataset' is a 3D+time sequence of volumes\n"
139 " ++ This must be a single imaging run -- that is, no discontinuities\n"
140 " in time from 3dTcat-ing multiple datasets together.\n"
141 "\n"
142 "* fbot = lowest frequency in the passband, in Hz\n"
143 " ++ fbot can be 0 if you want to do a lowpass filter only;\n"
144 " HOWEVER, the mean and Nyquist freq are always removed.\n"
145 "\n"
146 "* ftop = highest frequency in the passband (must be > fbot)\n"
147 " ++ if ftop > Nyquist freq, then it's a highpass filter only.\n"
148 "\n"
149 "* Set fbot=0 and ftop=99999 to do an 'allpass' filter.\n"
150 " ++ Except for removal of the 0 and Nyquist frequencies, that is.\n"
151 "\n"
152 "* You cannot construct a 'notch' filter with this program!\n"
153 " ++ You could use 3dRSFC followed by 3dcalc to get the same effect.\n"
154 " ++ If you are understand what you are doing, that is.\n"
155 " ++ Of course, that is the AFNI way -- if you don't want to\n"
156 " understand what you are doing, use Some other PrograM, and\n"
157 " you can still get Fine StatisticaL maps.\n"
158 "\n"
159 "* 3dRSFC will fail if fbot and ftop are too close for comfort.\n"
160 " ++ Which means closer than one frequency grid step df,\n"
161 " where df = 1 / (nfft * dt) [of course]\n"
162 "\n"
163 "* The actual FFT length used will be printed, and may be larger\n"
164 " than the input time series length for the sake of efficiency.\n"
165 " ++ The program will use a power-of-2, possibly multiplied by\n"
166 " a power of 3 and/or 5 (up to and including the 3rd power of\n"
167 " each of these: 3, 9, 27, and 5, 25, 125).\n"
168 "\n"
169 "* Note that the results of combining 3dDetrend and 3dRSFC will\n"
170 " depend on the order in which you run these programs. That's why\n"
171 " 3dRSFC has the '-ort' and '-dsort' options, so that the\n"
172 " time series filtering can be done properly, in one place.\n"
173 "\n"
174 "* The output dataset is stored in float format.\n"
175 "\n"
176 "* The order of processing steps is the following (most are optional), and\n"
177 " for the LFFs, the bandpass is done between the specified fbot and ftop,\n"
178 " while for the `whole spectrum' (i.e., fALFF denominator) the bandpass is:\n"
179 " done only to exclude the time series mean and the Nyquist frequency:\n"
180 " (0) Check time series for initial transients [does not alter data]\n"
181 " (1) Despiking of each time series\n"
182 " (2) Removal of a constant+linear+quadratic trend in each time series\n"
183 " (3) Bandpass of data time series\n"
184 " (4) Bandpass of -ort time series, then detrending of data\n"
185 " with respect to the -ort time series\n"
186 " (5) Bandpass and de-orting of the -dsort dataset,\n"
187 " then detrending of the data with respect to -dsort\n"
188 " (6) Blurring inside the mask [might be slow]\n"
189 " (7) Local PV calculation [WILL be slow!]\n"
190 " (8) L2 normalization [will be fast.]\n"
191 " (9) Calculate spectrum and amplitudes, for RSFC parameters.\n"
192 "\n"
193 "* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\n"
194 "--------\n"
195 "OPTIONS:\n"
196 "--------\n"
197 " -despike = Despike each time series before other processing.\n"
198 " ++ Hopefully, you don't actually need to do this,\n"
199 " which is why it is optional.\n"
200 " -ort f.1D = Also orthogonalize input to columns in f.1D\n"
201 " ++ Multiple '-ort' options are allowed.\n"
202 " -dsort fset = Orthogonalize each voxel to the corresponding\n"
203 " voxel time series in dataset 'fset', which must\n"
204 " have the same spatial and temporal grid structure\n"
205 " as the main input dataset.\n"
206 " ++ At present, only one '-dsort' option is allowed.\n"
207 " -nodetrend = Skip the quadratic detrending of the input that\n"
208 " occurs before the FFT-based bandpassing.\n"
209 " ++ You would only want to do this if the dataset\n"
210 " had been detrended already in some other program.\n"
211 " -dt dd = set time step to 'dd' sec [default=from dataset header]\n"
212 " -nfft N = set the FFT length to 'N' [must be a legal value]\n"
213 " -norm = Make all output time series have L2 norm = 1\n"
214 " ++ i.e., sum of squares = 1\n"
215 " -mask mset = Mask dataset\n"
216 " -automask = Create a mask from the input dataset\n"
217 " -blur fff = Blur (inside the mask only) with a filter\n"
218 " width (FWHM) of 'fff' millimeters.\n"
219 " -localPV rrr = Replace each vector by the local Principal Vector\n"
220 " (AKA first singular vector) from a neighborhood\n"
221 " of radius 'rrr' millimiters.\n"
222 " ++ Note that the PV time series is L2 normalized.\n"
223 " ++ This option is mostly for Bob Cox to have fun with.\n"
224 "\n"
225 " -input dataset = Alternative way to specify input dataset.\n"
226 " -band fbot ftop = Alternative way to specify passband frequencies.\n"
227 "\n"
228 " -prefix ppp = Set prefix name of output dataset. Name of filtered time\n"
229 " series would be, e.g., ppp_LFF+orig.*, and the parameter\n"
230 " outputs are named with obvious suffices.\n"
231 " -quiet = Turn off the fun and informative messages. (Why?)\n"
232 " -no_rs_out = Don't output processed time series-- just output\n"
233 " parameters (not recommended, since the point of\n"
234 " calculating RSFC params here is to have them be quite\n"
235 " related to the time series themselves which are used for\n"
236 " further analysis).\n"
237 " -un_bp_out = Output the un-bandpassed series as well (default is not \n"
238 " to). Name would be, e.g., ppp_unBP+orig.* .\n"
239 " with suffix `_unBP'.\n"
240 " -no_rsfa = If you don't want RSFA output (default is to do so).\n"
241 " -bp_at_end = A (probably unnecessary) switch to have bandpassing be \n"
242 " the very last processing step that is done in the\n"
243 " sequence of steps listed above; at Step 3 above, only \n"
244 " the time series mean and nyquist are BP'ed out, and then\n"
245 " the LFF series is created only after Step 9. NB: this \n"
246 " probably makes only very small changes for most\n"
247 " processing sequences (but maybe not, depending usage).\n"
248 "\n"
249 " -notrans = Don't check for initial positive transients in the data:\n"
250 " *OR* ++ The test is a little slow, so skipping it is OK,\n"
251 " -nosat if you KNOW the data time series are transient-free.\n"
252 " ++ Or set AFNI_SKIP_SATCHECK to YES.\n"
253 " ++ Initial transients won't be handled well by the\n"
254 " bandpassing algorithm, and in addition may seriously\n"
255 " contaminate any further processing, such as inter-\n"
256 " voxel correlations via InstaCorr.\n"
257 " ++ No other tests are made [yet] for non-stationary \n"
258 " behavior in the time series data.\n"
259 "\n"
260 "* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\n"
261 "\n"
262 " If you use this program, please reference the introductory/description\n"
263 " paper for the FATCAT toolbox:\n"
264 " Taylor PA, Saad ZS (2013). FATCAT: (An Efficient) Functional\n"
265 " And Tractographic Connectivity Analysis Toolbox. Brain \n"
266 " Connectivity 3(5):523-535.\n"
267 "____________________________________________________________________________\n"
268 );
269 PRINT_AFNI_OMP_USAGE(
270 " 3dRSFC" ,
271 " * At present, the only part of 3dRSFC that is parallelized is the\n"
272 " '-blur' option, which processes each sub-brick independently.\n"
273 ) ;
274 PRINT_COMPILE_DATE ; exit(0) ;
275 }
276
277 /*-- startup --*/
278
279 mainENTRY("3dRSFC"); machdep();
280 AFNI_logger("3dRSFC",argc,argv);
281 PRINT_VERSION("3dRSFC (from 3dBandpass by RW Cox): version THETA");
282 AUTHOR("PA Taylor");
283
284 nosat = AFNI_yesenv("AFNI_SKIP_SATCHECK") ;
285
286 nopt = 1 ;
287 while( nopt < argc && argv[nopt][0] == '-' ){
288
289 if( strcmp(argv[nopt],"-despike") == 0 ){ /* 08 Oct 2010 */
290 do_despike++ ; nopt++ ; continue ;
291 }
292
293 if( strcmp(argv[nopt],"-nfft") == 0 ){
294 int nnup ;
295 if( ++nopt >= argc ) ERROR_exit("need an argument after -nfft!") ;
296 nfft = (int)strtod(argv[nopt],NULL) ;
297 nnup = csfft_nextup_even(nfft) ;
298 if( nfft < 16 || nfft != nnup )
299 ERROR_exit("value %d after -nfft is illegal! Next legal value = %d",nfft,nnup) ;
300 nopt++ ; continue ;
301 }
302
303 if( strcmp(argv[nopt],"-blur") == 0 ){
304 if( ++nopt >= argc ) ERROR_exit("need an argument after -blur!") ;
305 blur = strtod(argv[nopt],NULL) ;
306 if( blur <= 0.0f ) WARNING_message("non-positive blur?!") ;
307 nopt++ ; continue ;
308 }
309
310 if( strcmp(argv[nopt],"-localPV") == 0 ){
311 if( ++nopt >= argc ) ERROR_exit("need an argument after -localpv!") ;
312 pvrad = strtod(argv[nopt],NULL) ;
313 if( pvrad <= 0.0f ) WARNING_message("non-positive -localpv?!") ;
314 nopt++ ; continue ;
315 }
316
317 if( strcmp(argv[nopt],"-prefix") == 0 ){
318 if( ++nopt >= argc ) ERROR_exit("need an argument after -prefix!") ;
319 prefix = strdup(argv[nopt]) ;
320 if( !THD_filename_ok(prefix) ) ERROR_exit("bad -prefix option!") ;
321 nopt++ ; continue ;
322 }
323
324 if( strcmp(argv[nopt],"-automask") == 0 ){
325 if( mask != NULL ) ERROR_exit("Can't use -mask AND -automask!") ;
326 do_automask = 1 ; nopt++ ; continue ;
327 }
328
329 if( strcmp(argv[nopt],"-mask") == 0 ){
330 THD_3dim_dataset *mset ;
331 if( ++nopt >= argc ) ERROR_exit("Need argument after '-mask'") ;
332 if( mask != NULL || do_automask ) ERROR_exit("Can't have two mask inputs") ;
333 mset = THD_open_dataset( argv[nopt] ) ;
334 CHECK_OPEN_ERROR(mset,argv[nopt]) ;
335 DSET_load(mset) ; CHECK_LOAD_ERROR(mset) ;
336 mask_nx = DSET_NX(mset); mask_ny = DSET_NY(mset); mask_nz = DSET_NZ(mset);
337 mask = THD_makemask( mset , 0 , 0.5f, 0.0f ) ; DSET_delete(mset) ;
338 if( mask == NULL ) ERROR_exit("Can't make mask from dataset '%s'",argv[nopt]) ;
339 nmask = THD_countmask( mask_nx*mask_ny*mask_nz , mask ) ;
340 if( verb ) INFO_message("Number of voxels in mask = %d",nmask) ;
341 if( nmask < 1 ) ERROR_exit("Mask is too small to process") ;
342 nopt++ ; continue ;
343 }
344
345 if( strcmp(argv[nopt],"-norm") == 0 ){
346 do_norm = 1 ; nopt++ ; continue ;
347 }
348
349 if( strcmp(argv[nopt],"-quiet") == 0 ){
350 verb = 0 ; nopt++ ; continue ;
351 }
352
353 if( strcmp(argv[nopt],"-no_rs_out") == 0 ){ // @@
354 SERIES_OUT = 0 ; nopt++ ; continue ;
355 }
356
357 if( strcmp(argv[nopt],"-un_bp_out") == 0 ){ // @@
358 UNBP_OUT = 1 ; nopt++ ; continue ;
359 }
360
361 if( strcmp(argv[nopt],"-no_rsfa") == 0 ){ // @@
362 DO_RSFA = 0 ; nopt++ ; continue ;
363 }
364
365 if( strcmp(argv[nopt],"-bp_at_end") == 0 ){ // @@
366 BP_LAST = 1 ; nopt++ ; continue ;
367 }
368
369
370
371
372 if( strcmp(argv[nopt],"-notrans") == 0 || strcmp(argv[nopt],"-nosat") == 0 ){
373 nosat = 1 ; nopt++ ; continue ;
374 }
375
376 if( strcmp(argv[nopt],"-ort") == 0 ){
377 if( ++nopt >= argc ) ERROR_exit("need an argument after -ort!") ;
378 if( ortar == NULL ) INIT_IMARR(ortar) ;
379 ortim = mri_read_1D( argv[nopt] ) ;
380 if( ortim == NULL ) ERROR_exit("can't read from -ort '%s'",argv[nopt]) ;
381 mri_add_name(argv[nopt],ortim) ;
382 ADDTO_IMARR(ortar,ortim) ;
383 nopt++ ; continue ;
384 }
385
386 if( strcmp(argv[nopt],"-dsort") == 0 ){
387 THD_3dim_dataset *qset ;
388 if( ++nopt >= argc ) ERROR_exit("need an argument after -dsort!") ;
389 if( nortset > 0 ) ERROR_exit("only 1 -dsort option is allowed!") ;
390 qset = THD_open_dataset(argv[nopt]) ;
391 CHECK_OPEN_ERROR(qset,argv[nopt]) ;
392 ortset = (THD_3dim_dataset **)realloc(ortset,
393 sizeof(THD_3dim_dataset *)*(nortset+1)) ;
394 ortset[nortset++] = qset ;
395 nopt++ ; continue ;
396 }
397
398 if( strncmp(argv[nopt],"-nodetrend",6) == 0 ){
399 qdet = 0 ; nopt++ ; continue ;
400 }
401
402 if( strcmp(argv[nopt],"-dt") == 0 ){
403 if( ++nopt >= argc ) ERROR_exit("need an argument after -dt!") ;
404 dt = (float)strtod(argv[nopt],NULL) ;
405 if( dt <= 0.0f ) WARNING_message("value after -dt illegal!") ;
406 nopt++ ; continue ;
407 }
408
409 if( strcmp(argv[nopt],"-input") == 0 ){
410 if( inset != NULL ) ERROR_exit("Can't have 2 -input options!") ;
411 if( ++nopt >= argc ) ERROR_exit("need an argument after -input!") ;
412 inset = THD_open_dataset(argv[nopt]) ;
413 CHECK_OPEN_ERROR(inset,argv[nopt]) ;
414
415 nopt++ ; continue ;
416 }
417
418 if( strncmp(argv[nopt],"-band",5) == 0 ){
419 if( ++nopt >= argc-1 ) ERROR_exit("need 2 arguments after -band!") ;
420 if( have_freq ) WARNING_message("second -band option replaces first one!") ;
421 fbot = strtod(argv[nopt++],NULL) ;
422 ftop = strtod(argv[nopt++],NULL) ;
423 have_freq = 1 ; continue ;
424 }
425
426 ERROR_exit("Unknown option: '%s'",argv[nopt]) ;
427 }
428
429 /** check inputs for reasonablositiness **/
430
431 if( !have_freq ){
432 if( nopt+1 >= argc )
433 ERROR_exit("Need frequencies on command line after options!") ;
434 fbot = (float)strtod(argv[nopt++],NULL) ;
435 ftop = (float)strtod(argv[nopt++],NULL) ;
436 }
437
438 if( inset == NULL ){
439 if( nopt >= argc )
440 ERROR_exit("Need input dataset name on command line after options!") ;
441 inset = THD_open_dataset(argv[nopt]) ;
442 CHECK_OPEN_ERROR(inset,argv[nopt]) ;
443
444 nopt++ ;
445 }
446 DSET_UNMSEC(inset) ;
447
448 if( fbot < 0.0f ) ERROR_exit("fbot value can't be negative!") ;
449 if( ftop <= fbot ) ERROR_exit("ftop value %g must be greater than fbot value %g!",ftop,fbot) ;
450
451 ntime = DSET_NVALS(inset) ;
452 if( ntime < 9 ) ERROR_exit("Input dataset is too short!") ;
453
454 if( nfft <= 0 ){
455 nfft = csfft_nextup_even(ntime) ;
456 if( verb ) INFO_message("Data length = %d FFT length = %d",ntime,nfft) ;
457 (void)THD_bandpass_set_nfft(nfft) ;
458 } else if( nfft < ntime ){
459 ERROR_exit("-nfft %d is less than data length = %d",nfft,ntime) ;
460 } else {
461 kk = THD_bandpass_set_nfft(nfft) ;
462 if( kk != nfft && verb )
463 INFO_message("Data length = %d FFT length = %d",ntime,kk) ;
464 }
465
466 if( dt <= 0.0f ){
467 dt = DSET_TR(inset) ;
468 if( dt <= 0.0f ){
469 WARNING_message("Setting dt=1.0 since input dataset lacks a time axis!") ;
470 dt = 1.0f ;
471 }
472 }
473 ftopALL = 1./dt ;// Aug,2016: should solve problem of a too-large
474 // value for THD_bandpass_vectors(), while still
475 // being >f_{Nyquist}
476
477 if( !THD_bandpass_OK(ntime,dt,fbot,ftop,1) ) ERROR_exit("Can't continue!") ;
478
479 nx = DSET_NX(inset); ny = DSET_NY(inset); nz = DSET_NZ(inset); nvox = nx*ny*nz;
480
481 /* check mask, or create it */
482
483 if( verb ) INFO_message("Loading input dataset time series" ) ;
484 DSET_load(inset) ;
485
486 if( mask != NULL ){
487 if( mask_nx != nx || mask_ny != ny || mask_nz != nz )
488 ERROR_exit("-mask dataset grid doesn't match input dataset") ;
489
490 } else if( do_automask ){
491 mask = THD_automask( inset ) ;
492 if( mask == NULL )
493 ERROR_message("Can't create -automask from input dataset?") ;
494 nmask = THD_countmask( DSET_NVOX(inset) , mask ) ;
495 if( verb ) INFO_message("Number of voxels in automask = %d",nmask);
496 if( nmask < 1 ) ERROR_exit("Automask is too small to process") ;
497
498 } else {
499 mask = (byte *)malloc(sizeof(byte)*nvox) ; nmask = nvox ;
500 memset(mask,1,sizeof(byte)*nvox) ;
501 // if( verb ) // @@ alert if aaaalllllll vox are going to be analyzed!
502 INFO_message("No mask ==> processing all %d voxels",nvox);
503 }
504
505 /* A simple check of dataset quality [08 Feb 2010] */
506
507 if( !nosat ){
508 float val ;
509 INFO_message(
510 "Checking dataset for initial transients [use '-notrans' to skip this test]") ;
511 val = THD_saturation_check(inset,mask,0,0) ; kk = (int)(val+0.54321f) ;
512 if( kk > 0 )
513 ININFO_message(
514 "Looks like there %s %d non-steady-state initial time point%s :-(" ,
515 ((kk==1) ? "is" : "are") , kk , ((kk==1) ? " " : "s") ) ;
516 else if( val > 0.3210f ) /* don't ask where this threshold comes from! */
517 ININFO_message(
518 "MAYBE there's an initial positive transient of 1 point, but it's hard to tell\n") ;
519 else
520 ININFO_message("No widespread initial positive transient detected :-)") ;
521 }
522
523 /* check -dsort inputs for match to inset */
524
525 for( kk=0 ; kk < nortset ; kk++ ){
526 if( DSET_NX(ortset[kk]) != nx ||
527 DSET_NY(ortset[kk]) != ny ||
528 DSET_NZ(ortset[kk]) != nz ||
529 DSET_NVALS(ortset[kk]) != ntime )
530 ERROR_exit("-dsort %s doesn't match input dataset grid" ,
531 DSET_BRIKNAME(ortset[kk]) ) ;
532 }
533
534 /* convert input dataset to a vectim, which is more fun */
535
536 // @@ convert BP'ing ftop/bot into indices for the DFT (below)
537 delf = 1.0/(ntime*dt);
538 ind_low = (int) rint(fbot/delf);
539 ind_high = (int) rint(ftop/delf);
540 if( ntime % 2 ) // nyquist number
541 N_ny = (ntime-1)/2;
542 else
543 N_ny = ntime/2;
544 sqnt = sqrt(ntime);
545 nt_fac = sqrt(ntime*(ntime-1));
546
547 // @@ if BP_LAST==0:
548 // now we go through twice, doing LFF bandpass for NumDen==0 and
549 // `full spectrum' processing for NumDen==1.
550 // if BP_LAST==1:
551 // now we go through once, doing only `full spectrum' processing
552 for( NumDen=0 ; NumDen<2 ; NumDen++) {
553 //if( NumDen==1 ){ // full spectrum
554 // fbot = fbotALL;
555 // ftop = ftopALL;
556 //}
557
558 // essentially, just doesn't BP here, and the perfect filtering at end
559 // is used for both still; this makes the final output spectrum
560 // contain only frequencies in range of 0.01-0.08
561 if( BP_LAST==1 )
562 INFO_message("Only doing filtering to LFFs at end!");
563
564
565 mrv = THD_dset_to_vectim( inset , mask , 0 ) ;
566 if( mrv == NULL ) ERROR_exit("Can't load time series data!?") ;
567 if( NumDen==1 )
568 DSET_unload(inset) ; // @@ only unload on 2nd pass
569
570 /* similarly for the ort vectors */
571
572 if( ortar != NULL ){
573 for( kk=0 ; kk < IMARR_COUNT(ortar) ; kk++ ){
574 ortim = IMARR_SUBIM(ortar,kk) ;
575 if( ortim->nx < ntime )
576 ERROR_exit("-ort file %s is shorter than input dataset time series",
577 ortim->name ) ;
578 ort = (float **)realloc( ort , sizeof(float *)*(nort+ortim->ny) ) ;
579 for( vv=0 ; vv < ortim->ny ; vv++ )
580 ort[nort++] = MRI_FLOAT_PTR(ortim) + ortim->nx * vv ;
581 }
582 }
583
584 /* all the real work now */
585
586 if( do_despike ){
587 int_pair nsp ;
588 if( verb ) INFO_message("Testing data time series for spikes") ;
589 nsp = THD_vectim_despike9( mrv ) ;
590 if( verb ) ININFO_message(" -- Squashed %d spikes from %d voxels",nsp.j,nsp.i) ;
591 }
592
593 if( verb ) INFO_message("Bandpassing data time series") ;
594
595 if( (BP_LAST==0) && (NumDen==0) )
596 (void)THD_bandpass_vectim( mrv , dt,fbot,ftop , qdet , nort,ort ) ;
597 else
598 (void)THD_bandpass_vectim( mrv , dt,fbotALL,ftopALL, qdet,nort,ort ) ;
599
600 /* OK, maybe a little more work */
601
602 if( nortset == 1 ){
603 MRI_vectim *orv ;
604 orv = THD_dset_to_vectim( ortset[0] , mask , 0 ) ;
605 if( orv == NULL ){
606 ERROR_message("Can't load -dsort %s",DSET_BRIKNAME(ortset[0])) ;
607 } else {
608 float *dp , *mvv , *ovv , ff ;
609 if( verb ) INFO_message("Orthogonalizing to bandpassed -dsort") ;
610 //(void)THD_bandpass_vectim( orv , dt,fbot,ftop , qdet , nort,ort ) ; //@@
611 if( (BP_LAST==0) && (NumDen==0) )
612 (void)THD_bandpass_vectim(orv,dt,fbot,ftop,qdet,nort,ort);
613 else
614 (void)THD_bandpass_vectim(orv,dt,fbotALL,ftopALL,qdet,nort,ort);
615
616 THD_vectim_normalize( orv ) ;
617 dp = malloc(sizeof(float)*mrv->nvec) ;
618 THD_vectim_vectim_dot( mrv , orv , dp ) ;
619 for( vv=0 ; vv < mrv->nvec ; vv++ ){
620 ff = dp[vv] ;
621 if( ff != 0.0f ){
622 mvv = VECTIM_PTR(mrv,vv) ; ovv = VECTIM_PTR(orv,vv) ;
623 for( kk=0 ; kk < ntime ; kk++ ) mvv[kk] -= ff*ovv[kk] ;
624 }
625 }
626 VECTIM_destroy(orv) ; free(dp) ;
627 }
628 }
629
630 if( blur > 0.0f ){
631 if( verb )
632 INFO_message("Blurring time series data spatially; FWHM=%.2f",blur) ;
633 mri_blur3D_vectim( mrv , blur ) ;
634 }
635 if( pvrad > 0.0f ){
636 if( verb )
637 INFO_message("Local PV-ing time series data spatially; radius=%.2f",pvrad) ;
638 THD_vectim_normalize( mrv ) ;
639 THD_vectim_localpv( mrv , pvrad ) ;
640 }
641 if( do_norm && pvrad <= 0.0f ){
642 if( verb ) INFO_message("L2 normalizing time series data") ;
643 THD_vectim_normalize( mrv ) ;
644 }
645
646 /* create output dataset, populate it, write it, then quit */
647 if( NumDen == 0 ) { // @@ BP'ed version; will do filt if BP_LAST
648
649 if(BP_LAST) // do bandpass here for BP_LAST
650 (void)THD_bandpass_vectim(mrv,dt,fbot,ftop,qdet,0,NULL);
651
652 if( verb ) INFO_message("Creating output dataset in memory, then writing it") ;
653 outset = EDIT_empty_copy(inset) ;
654 if(SERIES_OUT){
655 sprintf(out_lff,"%s_LFF",prefix);
656 EDIT_dset_items( outset , ADN_prefix,out_lff , ADN_none ) ;
657 tross_Copy_History( inset , outset ) ;
658 tross_Make_History( "3dBandpass" , argc,argv , outset ) ;
659 }
660 for( vv=0 ; vv < ntime ; vv++ )
661 EDIT_substitute_brick( outset , vv , MRI_float , NULL ) ;
662
663 #if 1
664 THD_vectim_to_dset( mrv , outset ) ;
665 #else
666 AFNI_OMP_START ;
667 #pragma omp parallel
668 { float *far , *var ; int *ivec=mrv->ivec ; int vv,kk ;
669 #pragma omp for
670 for( vv=0 ; vv < ntime ; vv++ ){
671 far = DSET_BRICK_ARRAY(outset,vv) ; var = mrv->fvec + vv ;
672 for( kk=0 ; kk < nmask ; kk++ ) far[ivec[kk]] = var[kk*ntime] ;
673 }
674 }
675 AFNI_OMP_END ;
676 #endif
677 VECTIM_destroy(mrv) ;
678 if(SERIES_OUT){ // @@
679 DSET_write(outset) ; if( verb ) WROTE_DSET(outset) ;
680 }
681 }
682 else{ // @@ non-BP'ed version
683 if( verb ) INFO_message("Creating output dataset 2 in memory") ;
684
685 // do this here because LFF version was also BP'ed at end.
686 if(BP_LAST) // do bandpass here for BP_LAST
687 (void)THD_bandpass_vectim(mrv,dt,fbotALL,ftopALL,qdet,0,NULL);
688
689 outsetALL = EDIT_empty_copy(inset) ;
690 if(UNBP_OUT){
691 sprintf(out_unBP,"%s_unBP",prefix);
692 EDIT_dset_items( outsetALL, ADN_prefix, out_unBP, ADN_none );
693 tross_Copy_History( inset , outsetALL ) ;
694 tross_Make_History( "3dRSFC" , argc,argv , outsetALL ) ;
695 }
696 for( vv=0 ; vv < ntime ; vv++ )
697 EDIT_substitute_brick( outsetALL , vv , MRI_float , NULL ) ;
698
699 #if 1
700 THD_vectim_to_dset( mrv , outsetALL ) ;
701 #else
702 AFNI_OMP_START ;
703 #pragma omp parallel
704 { float *far , *var ; int *ivec=mrv->ivec ; int vv,kk ;
705 #pragma omp for
706 for( vv=0 ; vv < ntime ; vv++ ){
707 far = DSET_BRICK_ARRAY(outsetALL,vv) ; var = mrv->fvec + vv ;
708 for( kk=0 ; kk < nmask ; kk++ ) far[ivec[kk]] = var[kk*ntime] ;
709 }
710 }
711 AFNI_OMP_END ;
712 #endif
713 VECTIM_destroy(mrv) ;
714 if(UNBP_OUT){
715 DSET_write(outsetALL) ; if( verb ) WROTE_DSET(outsetALL) ;
716 }
717 }
718 }// end of NumDen loop
719
720
721 // @@
722 INFO_message("Starting the (f)ALaFFel calcs") ;
723
724 // allocations
725 series1 = (double *)calloc(ntime,sizeof(double));
726 series2 = (double *)calloc(ntime,sizeof(double));
727 xx1 = (double *)calloc(2*ntime,sizeof(double));
728 xx2 = (double *)calloc(2*ntime,sizeof(double));
729 alff = (float *)calloc(nvox,sizeof(float));
730 malff = (float *)calloc(nvox,sizeof(float));
731 falff = (float *)calloc(nvox,sizeof(float));
732
733 if( (series1 == NULL) || (series2 == NULL)
734 || (xx1 == NULL) || (xx2 == NULL)
735 || (alff == NULL) || (malff == NULL) || (falff == NULL)) {
736 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
737 exit(122);
738 }
739 if(DO_RSFA) {
740 rsfa = (float *)calloc(nvox,sizeof(float));
741 mrsfa = (float *)calloc(nvox,sizeof(float));
742 frsfa = (float *)calloc(nvox,sizeof(float));
743 if( (rsfa == NULL) || (mrsfa == NULL) || (frsfa == NULL)) {
744 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
745 exit(123);
746 }
747 }
748
749
750 work = gsl_fft_real_workspace_alloc (ntime);
751 real1 = gsl_fft_real_wavetable_alloc (ntime);
752 real2 = gsl_fft_real_wavetable_alloc (ntime);
753 gsl_complex_packed_array compl_freqs1 = xx1;
754 gsl_complex_packed_array compl_freqs2 = xx2;
755
756
757
758
759 // *********************************************************************
760 // *********************************************************************
761 // ************** Falafelling = ALFF/fALFF calcs *****************
762 // *********************************************************************
763 // *********************************************************************
764
765 // Be now have the BP'ed data set (outset) and the non-BP'ed one
766 // (outsetALL). now we'll FFT both, get amplitudes in appropriate
767 // ranges, and calculate: ALFF, mALFF, fALFF,
768
769 ctr = 0;
770 for( kk=0; kk<nvox ; kk++) {
771 if(mask[kk]) {
772
773 // BP one, and unBP one, either for BP_LAST or !BP_LAST
774 for( m=0 ; m<ntime ; m++ ) {
775 series1[m] = THD_get_voxel(outset,kk,m);
776 series2[m] = THD_get_voxel(outsetALL,kk,m);
777 }
778
779
780 mm = gsl_fft_real_transform(series1, 1, ntime, real1, work);
781 mm = gsl_fft_halfcomplex_unpack(series1, compl_freqs1, 1, ntime);
782 mm = gsl_fft_real_transform(series2, 1, ntime, real2, work);
783 mm = gsl_fft_halfcomplex_unpack(series2, compl_freqs2, 1, ntime);
784
785 numer = 0.0f;
786 denom = 0.0f;
787 de_rsfa = 0.0f;
788 nu_rsfa = 0.0f;
789 for( m=1 ; m<N_ny ; m++ ) {
790 mm = 2*m;
791 pow2 = compl_freqs2[mm]*compl_freqs2[mm] +
792 compl_freqs2[mm+1]*compl_freqs2[mm+1]; // power
793 //pow2*=2;// factor of 2 since ampls are even funcs
794 denom+= (float) sqrt(pow2); // amplitude
795 de_rsfa+= (float) pow2;
796
797 if( ( m>=ind_low ) && ( m<=ind_high ) ){
798 pow1 = compl_freqs1[mm]*compl_freqs1[mm]+
799 compl_freqs1[mm+1]*compl_freqs1[mm+1];
800 //pow1*=2;
801 numer+= (float) sqrt(pow1);
802 nu_rsfa+= (float) pow1;
803 }
804 }
805
806 if( denom>0.000001 )
807 falff[kk] = numer/denom;
808 else
809 falff[kk] = 0.;
810 alff[kk] = 2*numer/sqnt;// factor of 2 since ampl is even funct
811 meanALFF+= alff[kk];
812
813 if(DO_RSFA){
814 nu_rsfa = sqrt(2*nu_rsfa); // factor of 2 since ampls
815 de_rsfa = sqrt(2*de_rsfa); // are even funcs
816 if( de_rsfa>0.000001 )
817 frsfa[kk] = nu_rsfa/de_rsfa;
818 else
819 frsfa[kk]=0.;
820 rsfa[kk] = nu_rsfa/nt_fac;
821 meanRSFA+= rsfa[kk];
822 }
823
824 ctr+=1;
825 }
826 }
827 meanALFF/= ctr;
828 meanRSFA/= ctr;
829
830 gsl_fft_real_wavetable_free(real1);
831 gsl_fft_real_wavetable_free(real2);
832 gsl_fft_real_workspace_free(work);
833
834 // ALFFs divided by mean of brain value
835 for( kk=0 ; kk<nvox ; kk++ )
836 if(mask[kk]){
837 malff[kk] = alff[kk]/meanALFF;
838 if(DO_RSFA)
839 mrsfa[kk] = rsfa[kk]/meanRSFA;
840 }
841 // **************************************************************
842 // **************************************************************
843 // Store and output
844 // **************************************************************
845 // **************************************************************
846
847 outsetALFF = EDIT_empty_copy( inset ) ;
848 sprintf(out_alff,"%s_ALFF",prefix);
849 EDIT_dset_items( outsetALFF,
850 ADN_nvals, 1,
851 ADN_datum_all , MRI_float ,
852 ADN_prefix , out_alff,
853 ADN_none ) ;
854 if( !THD_ok_overwrite() && THD_is_ondisk(DSET_HEADNAME(outsetALFF)) )
855 ERROR_exit("Can't overwrite existing dataset '%s'",
856 DSET_HEADNAME(outsetALFF));
857 EDIT_substitute_brick(outsetALFF, 0, MRI_float, alff);
858 alff=NULL;
859 THD_load_statistics(outsetALFF);
860 tross_Make_History("3dRSFC", argc, argv, outsetALFF);
861 THD_write_3dim_dataset(NULL, NULL, outsetALFF, True);
862
863 outsetfALFF = EDIT_empty_copy( inset ) ;
864 sprintf(out_falff,"%s_fALFF",prefix);
865 EDIT_dset_items( outsetfALFF,
866 ADN_nvals, 1,
867 ADN_datum_all , MRI_float ,
868 ADN_prefix , out_falff,
869 ADN_none ) ;
870 if( !THD_ok_overwrite() && THD_is_ondisk(DSET_HEADNAME(outsetfALFF)) )
871 ERROR_exit("Can't overwrite existing dataset '%s'",
872 DSET_HEADNAME(outsetfALFF));
873 EDIT_substitute_brick(outsetfALFF, 0, MRI_float, falff);
874 falff=NULL;
875 THD_load_statistics(outsetfALFF);
876 tross_Make_History("3dRSFC", argc, argv, outsetfALFF);
877 THD_write_3dim_dataset(NULL, NULL, outsetfALFF, True);
878
879
880
881 outsetmALFF = EDIT_empty_copy( inset ) ;
882 sprintf(out_malff,"%s_mALFF",prefix);
883 EDIT_dset_items( outsetmALFF,
884 ADN_nvals, 1,
885 ADN_datum_all , MRI_float ,
886 ADN_prefix , out_malff,
887 ADN_none ) ;
888 if( !THD_ok_overwrite() && THD_is_ondisk(DSET_HEADNAME(outsetmALFF)) )
889 ERROR_exit("Can't overwrite existing dataset '%s'",
890 DSET_HEADNAME(outsetmALFF));
891 EDIT_substitute_brick(outsetmALFF, 0, MRI_float, malff);
892 malff=NULL;
893 THD_load_statistics(outsetmALFF);
894 tross_Make_History("3dRSFC", argc, argv, outsetmALFF);
895 THD_write_3dim_dataset(NULL, NULL, outsetmALFF, True);
896
897 if(DO_RSFA){
898 outsetRSFA = EDIT_empty_copy( inset ) ;
899 sprintf(out_rsfa,"%s_RSFA",prefix);
900 EDIT_dset_items( outsetRSFA,
901 ADN_nvals, 1,
902 ADN_datum_all , MRI_float ,
903 ADN_prefix , out_rsfa,
904 ADN_none ) ;
905 if( !THD_ok_overwrite() && THD_is_ondisk(DSET_HEADNAME(outsetRSFA)) )
906 ERROR_exit("Can't overwrite existing dataset '%s'",
907 DSET_HEADNAME(outsetRSFA));
908 EDIT_substitute_brick(outsetRSFA, 0, MRI_float, rsfa);
909 rsfa=NULL;
910 THD_load_statistics(outsetRSFA);
911 tross_Make_History("3dRSFC", argc, argv, outsetRSFA);
912 THD_write_3dim_dataset(NULL, NULL, outsetRSFA, True);
913
914 outsetfRSFA = EDIT_empty_copy( inset ) ;
915 sprintf(out_frsfa,"%s_fRSFA",prefix);
916 EDIT_dset_items( outsetfRSFA,
917 ADN_nvals, 1,
918 ADN_datum_all , MRI_float ,
919 ADN_prefix , out_frsfa,
920 ADN_none ) ;
921 if( !THD_ok_overwrite() && THD_is_ondisk(DSET_HEADNAME(outsetfRSFA)) )
922 ERROR_exit("Can't overwrite existing dataset '%s'",
923 DSET_HEADNAME(outsetfRSFA));
924 EDIT_substitute_brick(outsetfRSFA, 0, MRI_float, frsfa);
925 frsfa=NULL;
926 THD_load_statistics(outsetfRSFA);
927 tross_Make_History("3dRSFC", argc, argv, outsetfRSFA);
928 THD_write_3dim_dataset(NULL, NULL, outsetfRSFA, True);
929
930 outsetmRSFA = EDIT_empty_copy( inset ) ;
931 sprintf(out_mrsfa,"%s_mRSFA",prefix);
932 EDIT_dset_items( outsetmRSFA,
933 ADN_nvals, 1,
934 ADN_datum_all , MRI_float ,
935 ADN_prefix , out_mrsfa,
936 ADN_none ) ;
937 if( !THD_ok_overwrite() && THD_is_ondisk(DSET_HEADNAME(outsetmRSFA)) )
938 ERROR_exit("Can't overwrite existing dataset '%s'",
939 DSET_HEADNAME(outsetmRSFA));
940 EDIT_substitute_brick(outsetmRSFA, 0, MRI_float, mrsfa);
941 mrsfa=NULL;
942 THD_load_statistics(outsetmRSFA);
943 tross_Make_History("3dRSFC", argc, argv, outsetmRSFA);
944 THD_write_3dim_dataset(NULL, NULL, outsetmRSFA, True);
945 }
946
947
948
949 // ************************************************************
950 // ************************************************************
951 // Freeing
952 // ************************************************************
953 // ************************************************************
954
955 DSET_delete(inset);
956 DSET_delete(outsetALL);
957 DSET_delete(outset);
958 DSET_delete(outsetALFF);
959 DSET_delete(outsetmALFF);
960 DSET_delete(outsetfALFF);
961 DSET_delete(outsetRSFA);
962 DSET_delete(outsetmRSFA);
963 DSET_delete(outsetfRSFA);
964
965 free(inset);
966 free(outsetALL);
967 free(outset);
968 free(outsetALFF);
969 free(outsetmALFF);
970 free(outsetfALFF);
971 free(outsetRSFA);
972 free(outsetmRSFA);
973 free(outsetfRSFA);
974
975 free(rsfa);
976 free(mrsfa);
977 free(frsfa);
978 free(alff);
979 free(malff);
980 free(falff);
981 free(mask);
982 free(series1);
983 free(series2);
984 free(xx1);
985 free(xx2);
986
987 exit(0) ;
988 }
989
990 /*----------------------------------------------------------------------------*/
991
THD_vectim_localpv(MRI_vectim * mrv,float rad)992 void THD_vectim_localpv( MRI_vectim *mrv , float rad )
993 {
994 unsigned short xran[3] = { 32701 , 22013 , 0x330e } ;
995 MCW_cluster *nbhd ;
996 int iv,kk,nn,nx,ny,nz,nxy , nind , *ind , xx,yy,zz , aa,bb,cc ;
997 float *pv , *fv ;
998 MRI_vectim *qrv ; /* workspace to hold results */
999
1000 nbhd = MCW_spheremask( mrv->dx,mrv->dy,mrv->dz , rad ) ;
1001 if( nbhd->num_pt <= 1 ){ THD_vectim_normalize(mrv); return; }
1002 ind = (int *)malloc(sizeof(int)*nbhd->num_pt) ;
1003 pv = (float *)malloc(sizeof(float)*mrv->nvals) ;
1004
1005 kk = thd_floatscan( mrv->nvec*mrv->nvals , mrv->fvec ) ;
1006 if( kk > 0 )
1007 WARNING_message("fixed %d float error%s before localPV",kk,(kk==1)?"\0":"s");
1008
1009 qrv = THD_vectim_copy(mrv) ; /* 08 Apr 2010: workspace */
1010
1011 nx = mrv->nx ; ny = mrv->ny ; nz = mrv->nz ; nxy = nx*ny ;
1012 for( iv=0 ; iv < mrv->nvec ; iv++ ){
1013 ind[0] = kk = mrv->ivec[iv] ; IJK_TO_THREE(kk,aa,bb,cc,nx,nxy) ;
1014 for( nind=nn=1 ; nn < nbhd->num_pt ; nn++ ){
1015 xx = aa + nbhd->i[nn] ; if( xx < 0 || xx >= nx ) continue ;
1016 yy = bb + nbhd->j[nn] ; if( yy < 0 || yy >= ny ) continue ;
1017 zz = cc + nbhd->k[nn] ; if( zz < 0 || zz >= nz ) continue ;
1018 ind[nind] = THREE_TO_IJK(xx,yy,zz,nx,nxy) ; nind++ ;
1019 }
1020
1021 nn = THD_vectim_subset_pv( mrv , nind,ind , pv , xran ) ;
1022 fv = VECTIM_PTR(qrv,iv) ; /* 08 Apr 2010: result goes in here, not mrv! */
1023 if( nn > 0 ){
1024 for( kk=0 ; kk < mrv->nvals ; kk++ ) fv[kk] = pv[kk] ;
1025 } else { /* should not happen */
1026 THD_normalize( mrv->nvals , fv ) ;
1027 }
1028 }
1029
1030 memcpy( mrv->fvec , qrv->fvec , sizeof(float)*mrv->nvec*mrv->nvals ) ;
1031 VECTIM_destroy(qrv) ;
1032
1033 kk = thd_floatscan( mrv->nvec*mrv->nvals , mrv->fvec ) ;
1034 if( kk > 0 )
1035 WARNING_message("fixed %d float error%s after localPV",kk,(kk==1)?"\0":"s");
1036
1037 KILL_CLUSTER(nbhd) ; free(ind) ; free(pv) ; return ;
1038 }
1039
1040 /*----------------------------------------------------------------------------*/
1041
THD_vectim_subset_pv(MRI_vectim * mrv,int nind,int * ind,float * ar,unsigned short xran[])1042 int THD_vectim_subset_pv( MRI_vectim *mrv, int nind, int *ind,
1043 float *ar, unsigned short xran[] )
1044 {
1045 int nvals , jj,kk,nkk ; register int ii ; float *fv , **xvec ;
1046
1047 if( mrv == NULL || nind <= 0 || ind == NULL || ar == NULL ) return 0 ;
1048
1049 nvals = mrv->nvals ;
1050 xvec = (float **)malloc(sizeof(float *)*nind) ;
1051
1052 for( nkk=jj=0 ; jj < nind ; jj++ ){
1053 kk = THD_vectim_ifind( ind[jj] , mrv ) ; if( kk < 0 ) continue ;
1054 xvec[nkk] = VECTIM_PTR(mrv,kk) ; nkk++ ;
1055 }
1056
1057 if( nkk == 0 ){ free(xvec) ; return 0 ; }
1058
1059 if( nkk == 1 ){
1060 for( ii=0 ; ii < nvals ; ii++ ) ar[ii] = xvec[0][ii] ;
1061 THD_normalize( nvals , ar ) ;
1062 free(xvec) ; return 1 ;
1063 }
1064
1065 (void)principal_vector( nvals,nkk , 1,xvec , ar , xvec[0] , NULL,xran ) ;
1066 return nkk ;
1067 }
1068