1 #include "mrilib.h"
2
3 /*----------------------------------------------------------------------------*/
4 /*-- This program implements and extends Anand Joshi's BrainSync algorithm. --*/
5 /*----------------------------------------------------------------------------*/
6
7 /*---------------------------*/
8 #undef MEMORY_CHECK
9 #if 0
10 #ifdef USING_MCW_MALLOC
11 # define MEMORY_CHECK(mm) \
12 do{ if( verb > 5 ) mcw_malloc_dump() ; \
13 if( verb > 1 ){ \
14 long long nb = mcw_malloc_total() ; \
15 if( nb > 0 ) INFO_message("Memory usage now = %s (%s): %s" , \
16 commaized_integer_string(nb) , \
17 approximate_number_string((double)nb) , (mm) ) ; \
18 ININFO_message(" status = %s",mcw_malloc_status(NULL,0)) ; \
19 } \
20 } while(0)
21 #endif
22 #endif
23
24 #ifndef MEMORY_CHECK
25 # define MEMORY_CHECK(mm) /*nada*/
26 #endif
27 /*---------------------------*/
28
29 /* various inputs from the options */
30 static int verb = 0 ;
31 static char *Qprefix = NULL , *Qbase = NULL ;
32 static char *Pprefix = NULL , *Pbase = NULL ;
33 static int do_norm = 0 ;
34
35 static int ct ; /* clock time in ms */
36 #define TIMER nice_time_string(NI_clock_time()-ct)+1
37
38 /* input datasets */
39 static THD_3dim_dataset *dsetB=NULL, *dsetC=NULL, *maskset=NULL ;
40
41 static int *bperm = NULL ; /* best permutation */
42 static THD_3dim_dataset *permset=NULL , *ortset=NULL ; /* output datasets */
43
44 #undef DO_DESPIKE /* at this time, despiking of inputs is disabled */
45 /* (expect that user will have done this already) */
46
47 /*----------------------------------------------------------------------------*/
48
BSY_help_the_pitiful_user(void)49 void BSY_help_the_pitiful_user(void) /* HELP ME IF YOU CAN, I'M FEELING DOWN */
50 {
51 printf(
52 "\n"
53 "Usage: 3dBrainSync [options]\n"
54 "\n"
55 "This program 'synchronizes' the -inset2 dataset to match the -inset1\n"
56 "dataset, as much as possible (average voxel-wise correlation), using the\n"
57 "same transformation on each input time series from -inset2:\n"
58 "\n"
59 " ++ With the -Qprefix option, the transformation is an orthogonal matrix,\n"
60 " computed as described in Joshi's original OHBM 2017 presentations,\n"
61 " and in the corresponding NeuroImage 2018 paper.\n"
62 " -->> Anand Joshi's presentation at OHBM was the genesis of this program.\n"
63 "\n"
64 " ++ With the -Pprefix option, the transformation is simply a\n"
65 " permutation of the time order of -inset2 (a very special case\n"
66 " of an orthogonal matrix).\n"
67 "\n"
68 " ++ The algorithms and a little discussion of the different features of\n"
69 " these two techniques are discussed in the METHODS section, infra.\n"
70 "\n"
71 " ++ At least one of '-Qprefix' or '-Pprefix' must be given, or\n"
72 " this program does not do anything! You can use both methods,\n"
73 " if you want to compare them.\n"
74 "\n"
75 " ++ 'Harmonize' might be a better name for what this program does,\n"
76 " but calling it 3dBrainHarm would probably not be good marketing\n"
77 " (except for Traumatic Brain Injury researchers?).\n"
78 "\n"
79 "One possible application of this program is to correlate resting state\n"
80 "FMRI datasets between subjects, voxel-by-voxel, as is sometimes done\n"
81 "with naturalistic stimuli (e.g., movie viewing).\n"
82 "\n"
83 "It would be amusing to see if within-subject resting state FMRI\n"
84 "runs can be BrainSync-ed better than between-subject runs.\n"
85 "\n"
86 "--------\n"
87 "OPTIONS:\n"
88 "--------\n"
89 " -inset1 dataset1 = Reference dataset\n"
90 "\n"
91 " -inset2 dataset2 = Dataset to be matched to the reference dataset,\n"
92 " as much as possible.\n"
93 " ++ These 2 datasets must be on the same spatial grid,\n"
94 " and must have the same number of time points!\n"
95 " ++ There must be at least twice as many voxels being\n"
96 " processed as there are time points (see '-mask', below).\n"
97 " ++ These are both MANDATORY 'options'.\n"
98 " ++ As usual in AFNI, since the computations herein are\n"
99 " voxel-wise, it is possible to input plain text .1D\n"
100 " files as datasets. When doing so, remember that\n"
101 " a ROW in the .1D file is interpreted as a time series\n"
102 " (single voxel's data). If your .1D files are oriented\n"
103 " so that time runs in down the COLUMNS, you will have to\n"
104 " transpose the inputs, which can be done on the command\n"
105 " line with the \\' operator, or externally using the\n"
106 " 1dtranspose program.\n"
107 " -->>++ These input datasets should be pre-processed first\n"
108 " to remove undesirable components (motions, baseline,\n"
109 " spikes, breathing, etc). Otherwise, you will be trying\n"
110 " to match artifacts between the datasets, which is not\n"
111 " likely to be interesting or useful. 3dTproject would be\n"
112 " one way to do this. Even better: afni_proc.py!\n"
113 " ++ In particular, the mean of each time series should have\n"
114 " been removed! Otherwise, the calculations are fairly\n"
115 " meaningless.\n"
116 "\n"
117 " -Qprefix qqq = Specifies the output dataset to be used for\n"
118 " the orthogonal matrix transformation.\n"
119 " ++ This will be the -inset2 dataset transformed\n"
120 " to be as correlated as possible (in time)\n"
121 " with the -inset1 dataset, given the constraint\n"
122 " that the transformation applied to each time\n"
123 " series is an orthogonal matrix.\n"
124 "\n"
125 " -Pprefix ppp = Specifies the output dataset to be used for\n"
126 " the permutation transformation.\n"
127 " ++ The output dataset is the -inset2 dataset\n"
128 " re-ordered in time, again to make the result\n"
129 " as correlated as possible with the -inset1\n"
130 " dataset.\n"
131 "\n"
132 " -normalize = Normalize the output dataset(s) so that each\n"
133 " time series has sum-of-squares = 1.\n"
134 " ++ This option is not usually needed in AFNI\n"
135 " (e.g., 3dTcorrelate does not care).\n"
136 "\n"
137 " -mask mset = Only operate on nonzero voxels in the mset dataset.\n"
138 " ++ Voxels outside the mask will not be used in computing\n"
139 " the transformation, but WILL be transformed for\n"
140 " your application and/or edification later.\n"
141 " ++ For FMRI purposes, a gray matter mask would make\n"
142 " sense here, or at least a brain mask.\n"
143 " ++ If no masking option is given, then all voxels\n"
144 " will be processed in computing the transformation.\n"
145 " This set will include all non-brain voxels (if any).\n"
146 " ++ Any voxel which is all constant in time\n"
147 " (in either input) will be removed from the mask.\n"
148 " ++ This mask dataset must be on the same spatial grid\n"
149 " as the other input datasets!\n"
150 "\n"
151 " -verb = Print some progress reports and auxiliary information.\n"
152 " ++ Use this option twice to get LOTS of progress\n"
153 " reports; mostly useful for debugging, or for fun.\n"
154 "\n"
155 "------\n"
156 "NOTES:\n"
157 "------\n"
158 "* Is this program useful? Not even The Shadow knows!\n"
159 " (But do NOT call it BS.)\n"
160 "\n"
161 "* The output dataset is in floating point format.\n"
162 "\n"
163 "* Although the goal of 3dBrainSync is to make the transformed\n"
164 " -inset2 as correlated (voxel-by-voxel) as possible with -inset1,\n"
165 " it does not actually compute or output that correlation dataset.\n"
166 " You can do that computation with program 3dTcorrelate, as in\n"
167 " 3dBrainSync -inset1 dataset1 -inset2 dataset2 \\\n"
168 " -Qprefix transformed-dataset2\n"
169 " 3dTcorrelate -polort -1 -prefix AB.pcor.nii \\\n"
170 " dataset1 transformed-dataset2\n"
171 "\n"
172 "* Besides the transformed dataset(s), if the '-verb' option is used,\n"
173 " some other (text formatted) files are written out:\n"
174 " {Qprefix}.sval.1D = singular values from the BC' decomposition\n"
175 " {Qprefix}.qmat.1D = Q matrix\n"
176 " {Pprefix}.perm.1D = permutation indexes p(i)\n"
177 " You probably do not have any use for these files; they are mostly\n"
178 " present to diagnose any problems.\n"
179 "\n"
180 "--------\n"
181 "METHODS:\n"
182 "--------\n"
183 "* Notation used in the explanations below:\n"
184 " M = Number of time points\n"
185 " N = Number of voxels > M (N = size of mask)\n"
186 " B = MxN matrix of time series from -inset1\n"
187 " C = MxN matrix of time series from -inset2\n"
188 " Both matrices will have each column normalized to\n"
189 " have sum-of-squares = 1 (L2 normalized) --\n"
190 " The program does this operation internally; you do not have\n"
191 " to ensure that the input datasets are so normalized.\n"
192 " Q = Desired orthgonal MxM matrix to transform C such that B-QC\n"
193 " is as small as possible (sum-of-squares = Frobenius norm).\n"
194 " That is, Q transforms dataset C to be as close as possible\n"
195 " to dataset B, given that Q is an orthogonal matrix.\n"
196 " normF(A) = sum_{ij} A_{ij}^2 = trace(AA') = trace(A'A).\n"
197 " NOTE: This norm is different from the matrix L2 norm.\n"
198 " NOTE: A' denotes the transpose of A.\n"
199 " NOTE: trace(A) = sum of diagonal element of square matrix A.\n"
200 " https://en.wikipedia.org/wiki/Matrix_norm\n"
201 "\n"
202 "* The expansion below shows why the matrix BC' is crucial to the analysis:\n"
203 " normF(B-QC) = trace( [B-QC][B'-C'Q'] )\n"
204 " = trace(BB') + trace(QCC'Q') - trace(BC'Q') - trace(QCB')\n"
205 " = trace(BB') + trace(C'C) - 2 trace(BC'Q')\n"
206 " The second term collapses because trace(AA') = trace(A'A), so\n"
207 " trace([QC][QC]') = trace([QC]'[QC]) = trace(C'Q'QC) = trace(C'C)\n"
208 " because Q is orthogonal. So the first 2 terms in the expansion of\n"
209 " normF(B-QC) do not depend on Q at all. Thus, to minimize normF(B-QC),\n"
210 " we have to maximize trace(BC'Q') = trace([B][QC]') = trace([QC][B]').\n"
211 "\n"
212 " Since the columns of B and C are the (normalized) time series,\n"
213 " each row represents the image at a particular time. So the (i,j)\n"
214 " element of BC' is the (spatial) dot product of the i-th TR image from\n"
215 " -inset1 with the j-th TR image from -inset2. Furthermore,\n"
216 " trace(BC') = trace(C'B) = sum of dot products (correlations)\n"
217 " of all time series. So maximizing trace(BC'Q') will maximize the\n"
218 " summed correlations of B (time series from -inset1) and QC\n"
219 " (transformed time series from -inset2).\n"
220 "\n"
221 " Note again that the sum of correlations (dot products) of all the time\n"
222 " series is equal to the sum of dot products of all the spatial images.\n"
223 " So the algorithm to find the transformation Q is to maximize the sum of\n"
224 " dot products of spatial images from B with Q-transformed spatial images\n"
225 " from C -- since there are fewer time points than voxels, this is more\n"
226 " efficient and elegant than trying to maximize the sum over voxels of dot\n"
227 " products of time series.\n"
228 "\n"
229 " If you use the '-verb' option, these summed correlations ('scores')\n"
230 " are printed to stderr during the analysis, for your fun and profit(?).\n"
231 "\n"
232 "*******************************************************************************\n"
233 "* Joshi method [-Qprefix]:\n"
234 " (a) compute MxM matrix B C'\n"
235 " (b) compute SVD of B C' = U S V' (U, S, V are MxM matrices)\n"
236 " (c) Q = U V'\n"
237 " [note: if B=C, then U=V, so Q=I, as it should]\n"
238 " (d) transform each time series from -inset2 using Q\n"
239 " This matrix Q is the solution to the restricted least squares\n"
240 " problem (i.e., restricted to have Q be an orthogonal matrix).\n"
241 " NOTE: The sum of the singular values in S is equal to the sum\n"
242 " of the time series dot products (correlations) in B and QC,\n"
243 " when Q is calculated as above.\n"
244 "\n"
245 " An article describing this method is available as:\n"
246 " AA Joshi, M Chong, RM Leahy.\n"
247 " Are you thinking what I'm thinking? Synchronization of resting fMRI\n"
248 " time-series across subjects.\n"
249 " NeuroImage v172:740-752 (2018).\n"
250 " https://doi.org/10.1016/j.neuroimage.2018.01.058\n"
251 " https://pubmed.ncbi.nlm.nih.gov/29428580/\n"
252 " https://www.google.com/search?q=joshi+brainsync\n"
253 "*******************************************************************************\n"
254 "* Permutation method [-Pprefix]:\n"
255 " (a) Compute B C' (same as above)\n"
256 " (b) Find a permutation p(i) of the integers {0..M-1} such\n"
257 " that sum_i { (BC')[i,p(i)] } is as large as possible\n"
258 " (i.e., p() is used as a permutation of the COLUMNS of BC').\n"
259 " This permutation is equivalent to post-multiplying BC'\n"
260 " by an orthogonal matrix P representing the permutation;\n"
261 " such a P is full of 0s except for a single 1 in each row\n"
262 " and each column.\n"
263 " (c) Permute the ROWS (time direction) of the time series matrix\n"
264 " from -inset2 using p().\n"
265 " Only an approximate (greedy) algorithm is used to find this\n"
266 " permutation; that is, the BEST permutation is not guaranteed to be found\n"
267 " (just a 'good' permutation -- it is the best thing I could code quickly :).\n"
268 "\n"
269 " Algorithm currently implemented (let D=BC' for notational simplicity):\n"
270 " 1) Find the largest element D(i,j) in the matrix.\n"
271 " Then the permutation at row i is p(i)=j.\n"
272 " Strike row i and column j out of the matrix D.\n"
273 " 2) Repeat, finding the largest element left, say at D(f,g).\n"
274 " Then p(f) = g. Strike row f and column g from the matrix.\n"
275 " Repeat until done.\n"
276 " (Choosing the largest possible element at each step is what makes this\n"
277 " method 'greedy'.) This permutation is not optimal but is pretty good,\n"
278 " and another step is used to improve it:\n"
279 " 3) For all pairs (i,j), p(i) and p(j) are swapped and that permutation\n"
280 " is tested to see if the trace gets bigger.\n"
281 " 4) This pair-wise swapping is repeated until it does not improve things\n"
282 " any more (typically, it improves the trace about 1-2%% -- not much).\n"
283 " The purpose of the pair swapping is to deal with situations where D looks\n"
284 " something like this: [ 1 70 ]\n"
285 " [ 70 99 ]\n"
286 " Step 1 would pick out 99, and Step 2 would pick out 1; that is,\n"
287 " p(2)=2 and then p(1)=1, for a total trace/score of 100. But swapping\n"
288 " 1 and 2 would give a total trace/score of 140. In practice, extreme versions\n"
289 " of this situation do not seem common with real FMRI data, probably because\n"
290 " the subject's brain isn't actively conspiring against this algorithm :)\n"
291 "\n"
292 " [Something called the 'Hungarian algorithm' can solve for the optimal]\n"
293 " [permutation exactly, but I've not had the inclination to program it.]\n"
294 "\n"
295 " This whole permutation optimization procedure is very fast: about 1 second.\n"
296 " In the RS-FMRI data I've tried this on, the average time series correlation\n"
297 " resulting from this optimization is 30-60%% of that which comes from\n"
298 " optimizing over ALL orthogonal matrices (Joshi method). If you use '-verb',\n"
299 " the stderr output line that looks like this\n"
300 " + corr scores: original=-722.5 Q matrix=22366.0 permutation=12918.7 57.8%%\n"
301 " shows trace(BC') before any transforms, with the Q matrix transform,\n"
302 " and with the permutation transform. As explained above, trace(BC') is\n"
303 " the summed correlations of the time series (since the columns of B and C\n"
304 " are normalized prior to the optimizations); in this example, the ratio of\n"
305 " the average time series correlation between the permutation method and the\n"
306 " Joshi method is about 58%% (in a gray matter mask with 72221 voxels).\n"
307 "\n"
308 "* Results from the permutation method MUST be less correlated (on average)\n"
309 " with -inset1 than the Joshi method's results: the permutation can be\n"
310 " thought of as an orthogonal matrix containing only 1s and 0s, and the BEST\n"
311 " possible orthogonal matrix, from Joshi's method, has more general entries.\n"
312 " ++ However, the permutation method has an obvious interpretation\n"
313 " (re-ordering time points), while the general method linearly combines\n"
314 " different time points (perhaps far apart); the interpretation of this\n"
315 " combination in terms of synchronizing brain activity is harder to intuit\n"
316 " (at least for me).\n"
317 " ++ Another feature of a permutation-only transformation is that it cannot\n"
318 " change the sign of data, unlike a general orthgonal matrix; e.g.,\n"
319 " [ 0 -1 0 ]\n"
320 " [-1 0 0 ]\n"
321 " [ 0 0 1 ], which swaps the first 2 time points AND negates them,\n"
322 " and leave the 3rd time point unchanged, is a valid orthogonal\n"
323 " matrix. For rs-FMRI datasets, this consideration might not be important,\n"
324 " since rs-FMRI correlations are generally positive, so don't often need\n"
325 " sign-flipping to make them so.\n"
326 "*******************************************************************************\n"
327 "\n"
328 "* This program is NOT multi-threaded. Typically, I/O is the biggest part of\n"
329 " the run time (at least, for the cases I've tested). The '-verb' option\n"
330 " will give progress reports with elapsed-time stamps, making it easy to\n"
331 " see which parts of the program take the most time.\n"
332 "\n"
333 "* Author: RWCox, servant of the ChronoSynclastic Infundibulum - July 2017\n"
334 "\n"
335 "* Thanks go to Anand Joshi for his clear exposition of BrainSync at OHBM 2017,\n"
336 " and his encouragement about the development of this program.\n"
337 ) ;
338
339 PRINT_COMPILE_DATE ; exit(0) ;
340 }
341
342 /*----------------------------------------------------------------------------*/
343
344 #undef A
345 #undef QQ
346 #define A(i,j) amat[(i)+(j)*m]
347 #define QQ(i,j) qmat[(i)+(j)*m]
348
349 #undef BIGG
350 #define BIGG 1.111e37f
351
352 /*-------- Find the "best" permutation, in
353 the sense of the largest trace of the row-permuted amat
354 (or at least, an approximation of the "best" permutation :)
355 Matrix qmat is used only for comparison, and can be NULL.
356 Input matrices qmat and amat are square, m x m. ---------*/
357
find_best_permutation(int m,float * qmat,double * amat)358 static int * find_best_permutation( int m , float *qmat , double *amat )
359 {
360 int nrand ;
361 float *rvec,*qvec , rcost=0.0f, qcost=0.0f , bestcost, rval ;
362 int *perm , *bestperm=NULL , ii,jj,rr,qq ;
363
364 if( m < 2 || amat == NULL ) return NULL ; /* should not happen */
365
366 ENTRY("find_best_permutation") ;
367
368 /* find score of the Q matrix, if it was input [for comparison] */
369
370 if( qmat != NULL ){
371 for( qcost=0.0f,jj=0 ; jj < m ; jj++ ){
372 for( ii=0 ; ii < m ; ii++ ) qcost += A(ii,jj)*QQ(ii,jj) ;
373 }
374 if( verb > 1 )
375 ININFO_message("-- permutation search: qmat score = %g [%s]",qcost,TIMER) ;
376 } else {
377 if( verb > 1 )
378 ININFO_message("-- permutation search starts [%s]",TIMER) ;
379 }
380
381 /* space for the trial permutation */
382
383 perm = (int *)malloc(sizeof(int)*m) ;
384 bestcost = -BIGG ; /* not very good */
385
386 /* try a very simple greedy algorithm */
387 #undef AA
388 #define AA(i,j) aamat[(i)+(j)*m]
389 { int kk,pi,pj , mmm=m*m , *jdone ;
390 float *aamat=(float *)malloc(sizeof(float)*m*m) , bbcost=0.0f ;
391
392 /* copy amat into aamat */
393 for( kk=0 ; kk < mmm ; kk++ ) aamat[kk] = (float)amat[kk] ;
394
395 bestperm = (int *)malloc(sizeof(int)*m) ; /* best yet */
396 jdone = (int *)malloc(sizeof(int)*m) ; /* cols already done */
397 for( ii=0 ; ii < m ; ii++ ){ bestperm[ii] = -1; jdone[ii] = 0; }
398
399 for( rr=0 ; rr < m ; rr++ ){ /* pick biggest element left */
400 rval = -BIGG ; pi = pj = -1 ;
401 for( jj=0 ; jj < m ; jj++ ){
402 if( jdone[jj] ) continue ; /* already did this column */
403 for( ii=0 ; ii < m ; ii++ ){
404 if( AA(ii,jj) > rval ){ rval = AA(ii,jj); pi=ii; pj=jj; }
405 }
406 }
407 if( pi < 0 ) break ; /* should not be possible */
408 bestperm[pi] = pj ; jdone[pj] = 1 ;
409 for( kk=0 ; kk < m ; kk++ ){ /* strike out this row+col */
410 AA(pi,kk) = AA(kk,pj) = -BIGG;
411 }
412 }
413 free(aamat); free(jdone);
414
415 /* check if every entry in bestperm was set above */
416 for( ii=0 ; ii < m && bestperm[ii] >= 0 ; ii++ ) ; /*nada*/
417 if( ii < m ){
418 free(bestperm); bestperm = NULL; /* should never happen */
419 } else { /* compute cost of bestperm */
420 for( rcost=0.0f,ii=0 ; ii < m ; ii++ ) rcost += A(ii,bestperm[ii]) ;
421 bestcost = bbcost = rcost ;
422 if( verb > 1 )
423 ININFO_message(" initial greedy perm score = %g",bbcost) ;
424 if( verb > 2 ){ /* debugggin output */
425 fprintf(stderr," perm:") ;
426 for( ii=0 ; ii < m ; ii++ )
427 fprintf(stderr," A[%d,%d]=%g",ii,bestperm[ii],A(ii,bestperm[ii])) ;
428 fprintf(stderr,"\n") ;
429 }
430 }
431
432 /* bestperm should not be NULL here, unless the Gods are out to get me */
433
434 if( bestperm != NULL ){ /* swap pairs, look for improvement */
435 int ntry=0, ndone=0, ntot=0 ;
436 do{ /* loopback to try another pairswap fun */
437 for( ndone=ii=0 ; ii < m-1 ; ii++ ){
438 for( jj=ii+1 ; jj < m ; jj++ ){ /* swapping ii,jj */
439 memcpy(perm,bestperm,sizeof(int)*m) ;
440 kk = perm[ii] ; perm[ii] = perm[jj] ; perm[jj] = kk ;
441 for( rcost=0.0f,kk=0 ; kk < m ; kk++ ) rcost += A(kk,perm[kk]) ;
442 if( rcost > bestcost ){ /* this is better? keep it! */
443 free(bestperm); bestperm = perm; bestcost = rcost;
444 perm = (int *)malloc(sizeof(int)*m); ndone++ ; ntot++ ;
445 if( verb > 1 )
446 ININFO_message(" improved perm score = %g swap [%d,%d]",rcost,ii,jj) ;
447 }
448 }} /* end of loops over ii,jj */
449 } while( ndone > 0 && ++ntry < 29 ) ; /* iterate, but don't try forever! */
450 if( verb > 1 && ntot > 0 ){ /* report on iteration results */
451 if( bbcost > 0.0f ) bbcost = 100.0f*(bestcost-bbcost)/bbcost ;
452 else bbcost = 0.0f ;
453 ININFO_message("-- finish permutation after %d swaps in %d iterates: %.2f%% improvement",
454 ntot,ntry,bbcost) ;
455 }
456 if( verb > 2 ){ /* debugging output */
457 fprintf(stderr," perm:") ;
458 for( ii=0 ; ii < m ; ii++ )
459 fprintf(stderr," A[%d,%d]=%g",ii,bestperm[ii],A(ii,bestperm[ii])) ;
460 fprintf(stderr,"\n") ;
461 }
462 }
463 }
464 #undef AA
465
466 if( bestperm == NULL ){ /* should never happen: use old method */
467 nrand = 9999*m ; /* of random trials (not good) */
468 rvec = (float *)malloc(sizeof(float)*m) ;
469 qvec = (float *)malloc(sizeof(float)*m) ;
470 WARNING_message("Greedy algorithm for permutation failed; using backup") ;
471 for( rr=0 ; rr < nrand ; rr++ ){
472 for( ii=0 ; ii < m ; ii++ ){
473 perm[ii] = ii; rvec[ii] = 0.0f;
474 qvec[ii] = (float)(drand48()+drand48()-1.0);
475 }
476 qsort_float(m,qvec) ;
477 for( jj=0 ; jj < m ; jj++ ){
478 rval = qvec[jj] ;
479 for( ii=0 ; ii < m ; ii++ ) rvec[ii] += QQ(ii,jj)*rval ;
480 }
481 qsort_floatint(m,rvec,perm) ;
482
483 for( rcost=0.0f,ii=0 ; ii < m ; ii++ ) rcost += A(ii,perm[ii]) ;
484 if( rcost > bestcost ){
485 if( bestperm != NULL ) free(bestperm) ;
486 bestperm = perm; bestcost = rcost;
487 perm = (int *)malloc(sizeof(int)*m);
488 if( verb > 1 )
489 ININFO_message(" new perm score = %g [%d]",rcost,rr) ;
490 }
491 }
492 free(rvec); free(qvec);
493 }
494
495 free(perm); /* only bestperm is our friend now */
496
497 if( verb && qcost > 0.0f ){
498 for( rcost=0.0f,ii=0 ; ii < m ; ii++ ) rcost += A(ii,ii) ;
499 ININFO_message("corr scores: "
500 "original=%.1f Q matrix=%.1f permutation=%.1f %.1f%%",
501 rcost, qcost, bestcost , 100.0f*bestcost/qcost ) ;
502 }
503
504 RETURN(bestperm) ;
505 }
506
507 #undef A
508 #undef QQ
509
510 /*----------------------------------------------------------------------------*/
511 /* Given m X n matrices bmat and cmat, compute the orthogonal m X m matrix
512 qmat that 'best' transforms cmat to bmat.
513 1) normalize all columns of bmat and cmat
514 2) compute m X m matrix amat = [bmat] * [cmat]'
515 3) SVD that to get amat = [umat] * [sigma] * [vmat]'
516 4) qmat = [umat] * [vmat]'
517 Also compute the best permutation while we're at it (that's fast).
518 *//*--------------------------------------------------------------------------*/
519
compute_joshi_matrix(int m,int n,float * bmat,float * cmat,float * qmat)520 static void compute_joshi_matrix( int m , int n ,
521 float *bmat , float *cmat , float *qmat )
522 {
523 int ii,jj,kk , kbot ;
524 double *bmatn , *cmatn ;
525 double *amat , *umat , *vmat , *sval ;
526 register double bsum , csum ;
527
528 ENTRY("compute_joshi_matrix") ;
529
530 /* temp matrices */
531
532 bmatn = (double *)calloc( sizeof(double) , m*n ) ; /* normalized bmat */
533 cmatn = (double *)calloc( sizeof(double) , m*n ) ; /* normalized cmat */
534 amat = (double *)calloc( sizeof(double) , m*m ) ; /* [bmatn] * [cmatn]' */
535
536 /* macros for 2D array indexing */
537
538 #undef B
539 #undef C
540 #undef A
541 #undef U
542 #undef V
543 #undef BN
544 #undef CN
545
546 #define B(i,j) bmat[(i)+(j)*m]
547 #define C(i,j) cmat[(i)+(j)*m]
548 #define BN(i,j) bmatn[(i)+(j)*m]
549 #define CN(i,j) cmatn[(i)+(j)*m]
550 #define A(i,j) amat[(i)+(j)*m]
551 #define U(i,j) umat[(i)+(j)*m]
552 #define V(i,j) vmat[(i)+(j)*m]
553 #define QQ(i,j) qmat[(i)+(j)*m]
554
555 /* copy input matrices into bmatn and cmatn, normalizing as we go */
556
557 if( verb > 1 ) ININFO_message("Normalizing time series [%s]",TIMER) ;
558
559 for( jj=0 ; jj < n ; jj++ ){
560 bsum = csum = 0.0 ;
561 for( ii=0 ; ii < m ; ii++ ){
562 bsum += B(ii,jj)*B(ii,jj) ;
563 csum += C(ii,jj)*C(ii,jj) ;
564 }
565 if( bsum > 0.0 ) bsum = 1.0 / sqrt(bsum) ; /* L2 normalizing */
566 if( csum > 0.0 ) csum = 1.0 / sqrt(csum) ; /* factors */
567 for( ii=0 ; ii < m ; ii++ ){
568 BN(ii,jj) = B(ii,jj)*bsum ;
569 CN(ii,jj) = C(ii,jj)*csum ;
570 }
571 }
572 MEMORY_CHECK("a") ;
573
574 /* form A matrix = BN * CN' */
575
576 if( verb > 1 ) ININFO_message("forming BC' matrix [%s]",TIMER) ;
577
578 #if 0 /* old method for computing BC' */
579 kbot = n%2 ; /* 1 if odd, 0 if even */
580 for( jj=0 ; jj < m ; jj++ ){
581 for( ii=0 ; ii < m ; ii++ ){
582 bsum = (kbot) ? BN(ii,0)*CN(jj,0) : 0.0 ;
583 for( kk=kbot ; kk < n ; kk+=2 ) /* unrolled by 2 */
584 bsum += BN(ii,kk)*CN(jj,kk) + BN(ii,kk+1)*CN(jj,kk+1) ;
585 A(ii,jj) = bsum ;
586 }}
587 #else /* new method (faster) */
588 for( jj=0 ; jj < m ; jj++ ){
589 for( ii=0 ; ii < m ; ii++ ) A(ii,jj) = 0.0 ;
590 }
591 for( kk=0 ; kk < n ; kk++ ){
592 for( jj=0 ; jj < m ; jj++ ){
593 bsum = CN(jj,kk) ;
594 for( ii=0 ; ii < m ; ii++ ) A(ii,jj) += BN(ii,kk)*bsum ;
595 }
596 }
597 #endif
598 MEMORY_CHECK("b") ;
599
600 free(bmatn) ; free(cmatn) ;
601 umat = (double *)calloc( sizeof(double),m*m ) ; /* SVD outputs */
602 vmat = (double *)calloc( sizeof(double),m*m ) ;
603 sval = (double *)calloc( sizeof(double),m ) ;
604
605 /* compute SVD of scaled matrix */
606
607 if( verb > 1 ) ININFO_message("SVD-ing BC' matrix %dx%d [%s]",m,m,TIMER) ;
608
609 svd_double( m , m , amat , sval , umat , vmat ) ;
610
611 /* since sval is otherwise unused, sort it so largest are firstest */
612
613 for( ii=0 ; ii < m ; ii++ ) sval[ii] = -sval[ii] ;
614 qsort_double( m , sval ) ;
615 for( ii=0 ; ii < m ; ii++ ) sval[ii] = -sval[ii] ;
616
617 bsum = 0.0 ;
618 for( ii=0 ; ii < m ; ii++ ) bsum += sval[ii] ;
619 if( verb > 1 )
620 ININFO_message("first 5 singular values: "
621 " %.2f [%.2f%%] "
622 " %.2f [%.2f%%] "
623 " %.2f [%.2f%%] "
624 " %.2f [%.2f%%] "
625 " %.2f [%.2f%%] ... sum = %.1f" ,
626 sval[0], 100.0*sval[0]/bsum ,
627 sval[1], 100.0*sval[1]/bsum ,
628 sval[2], 100.0*sval[2]/bsum ,
629 sval[3], 100.0*sval[3]/bsum ,
630 sval[4], 100.0*sval[4]/bsum , bsum ) ;
631
632 /* write the singular values to a file? [30 Jul 2017] */
633
634 if( Qbase != NULL && verb ){
635 char *fname ; MRI_IMAGE *qim ; float *qar ;
636 fname = (char *)malloc(sizeof(char)*(strlen(Qbase)+32)) ;
637 strcpy(fname,Qbase) ;
638 strcat(fname,".sval.1D") ;
639 qim = mri_new( m,1, MRI_float ) ; qar = MRI_FLOAT_PTR(qim) ;
640 for( ii=0 ; ii < m ; ii++ ) qar[ii] = sval[ii] ;
641 mri_write_1D( fname , qim ) ;
642 ININFO_message("wrote singular values to %s",fname) ;
643 mri_free(qim) ; free(fname) ;
644 }
645
646 /* compute QQ = output matrix = U V' */
647
648 if( verb > 1 ) ININFO_message("computing Q matrix [%s]",TIMER) ;
649
650 for( jj=0 ; jj < m ; jj++ ){
651 for( ii=0 ; ii < m ; ii++ ){
652 bsum = 0.0 ;
653 for( kk=0 ; kk < m ; kk++ ) bsum += U(ii,kk)*V(jj,kk) ;
654 QQ(ii,jj) = (float)bsum ;
655 }}
656 MEMORY_CHECK("d") ;
657
658 free(umat) ; free(vmat) ;
659
660 /* also find permutation to approximate qmat */
661
662 bperm = find_best_permutation( m , qmat , amat ) ;
663
664 free(amat) ; free(sval) ;
665
666 /* write the Q matrix to a file? */
667
668 if( Qbase != NULL && verb ){
669 char *fname ; MRI_IMAGE *qim ; float *qar ;
670 fname = (char *)malloc(sizeof(char)*(strlen(Qbase)+32)) ;
671 strcpy(fname,Qbase) ;
672 strcat(fname,".qmat.1D") ;
673 qim = mri_new_vol_empty( m,m,1, MRI_float ) ;
674 mri_fix_data_pointer( qmat , qim ) ;
675 mri_write_1D( fname , qim ) ;
676 ININFO_message("wrote Q matrix to %s",fname) ;
677 mri_clear_data_pointer(qim) ; mri_free(qim) ; free(fname) ;
678 }
679
680 /* write the permutation to a file? */
681
682 if( Pbase != NULL && bperm != NULL && verb ){
683 char *fname ; MRI_IMAGE *qim ; float *qar ;
684 fname = (char *)malloc(sizeof(char)*(strlen(Pbase)+32)) ;
685 strcpy(fname,Pbase) ;
686 strcat(fname,".perm.1D") ;
687 qim = mri_new(m,1,MRI_float) ; qar = MRI_FLOAT_PTR(qim) ;
688 for( ii=0 ; ii < m ; ii++ ) qar[ii] = bperm[ii] ;
689 mri_write_1D(fname,qim); mri_free(qim);
690 ININFO_message("wrote permutation vector to %s",fname) ;
691 free(fname) ;
692 }
693
694 EXRETURN ;
695 }
696
697 /*----------------------------------------------------------------------------*/
698
is_vector_constant(int n,float * v)699 static INLINE int is_vector_constant( int n , float *v )
700 {
701 int ii ;
702 for( ii=1 ; ii < n && v[ii] == v[0] ; ii++ ) ; /*nada*/
703 return (ii==n) ;
704 }
705
706 #ifdef DO_DESPIKE
707 /*------------------- 08 Oct 2010: functions for despiking ----------------*/
708
709 #undef SWAP
710 #define SWAP(x,y) (temp=x,x=y,y=temp)
711
712 #undef SORT2
713 #define SORT2(a,b) if(a>b) SWAP(a,b)
714
715 /*--- fast median of 9 values ---*/
716
median9f(float * p)717 static INLINE float median9f(float *p)
718 {
719 register float temp ;
720 SORT2(p[1],p[2]) ; SORT2(p[4],p[5]) ; SORT2(p[7],p[8]) ;
721 SORT2(p[0],p[1]) ; SORT2(p[3],p[4]) ; SORT2(p[6],p[7]) ;
722 SORT2(p[1],p[2]) ; SORT2(p[4],p[5]) ; SORT2(p[7],p[8]) ;
723 SORT2(p[0],p[3]) ; SORT2(p[5],p[8]) ; SORT2(p[4],p[7]) ;
724 SORT2(p[3],p[6]) ; SORT2(p[1],p[4]) ; SORT2(p[2],p[5]) ;
725 SORT2(p[4],p[7]) ; SORT2(p[4],p[2]) ; SORT2(p[6],p[4]) ;
726 SORT2(p[4],p[2]) ; return(p[4]) ;
727 }
728 #undef SORT2
729 #undef SWAP
730
731 /*--- get the local median and MAD of values vec[j-4 .. j+4] ---*/
732
733 #undef mead9
734 #define mead9(j) \
735 { float qqq[9] ; int jj = (j)-4 ; \
736 if( jj < 0 ) jj = 0; else if( jj+8 >= num ) jj = num-9; \
737 qqq[0] = vec[jj+0]; qqq[1] = vec[jj+1]; qqq[2] = vec[jj+2]; \
738 qqq[3] = vec[jj+3]; qqq[4] = vec[jj+4]; qqq[5] = vec[jj+5]; \
739 qqq[6] = vec[jj+6]; qqq[7] = vec[jj+7]; qqq[8] = vec[jj+8]; \
740 med = median9f(qqq); qqq[0] = fabsf(qqq[0]-med); \
741 qqq[1] = fabsf(qqq[1]-med); qqq[2] = fabsf(qqq[2]-med); \
742 qqq[3] = fabsf(qqq[3]-med); qqq[4] = fabsf(qqq[4]-med); \
743 qqq[5] = fabsf(qqq[5]-med); qqq[6] = fabsf(qqq[6]-med); \
744 qqq[7] = fabsf(qqq[7]-med); qqq[8] = fabsf(qqq[8]-med); \
745 mad = median9f(qqq); }
746
747 /*-------------------------------------------------------------------------*/
748 /*! Remove spikes from a time series, in a very simplistic way.
749 Return value is the number of spikes that were squashed [RWCox].
750 *//*-----------------------------------------------------------------------*/
751
despike9(int num,float * vec)752 static int despike9( int num , float *vec )
753 {
754 int ii , nsp ; float *zma,*zme , med,mad,val ;
755
756 if( num < 9 || vec == NULL ) return 0 ;
757 zme = (float *)malloc(sizeof(float)*num) ;
758 zma = (float *)malloc(sizeof(float)*num) ;
759
760 for( ii=0 ; ii < num ; ii++ ){
761 mead9(ii) ; zme[ii] = med ; zma[ii] = mad ;
762 }
763 mad = qmed_float(num,zma) ; free(zma) ;
764 if( mad <= 0.0f ){ free(zme); return 0; } /* should not happen */
765 mad *= 6.789f ; /* threshold value */
766
767 for( nsp=ii=0 ; ii < num ; ii++ )
768 if( fabsf(vec[ii]-zme[ii]) > mad ){ vec[ii] = zme[ii]; nsp++; }
769
770 free(zme) ; return nsp ;
771 }
772 #undef mead9
773 /*----------------------------------------------------------------------------*/
774 #endif /* DO_DESPIKE */
775
776 /*----------------------------------------------------------------------------*/
777 /* compute the output datasets */
778
BSY_process_data(void)779 void BSY_process_data(void)
780 {
781 byte *vmask ;
782 int nvmask , ii,jj,kk , mm , nvox , ncut=0 ;
783 MRI_IMAGE *bim, *cim ;
784 float *bar, *car, *bmat, *cmat, *qmat, *cvec ;
785 register float csum ;
786
787 ENTRY("BSY_process_data") ;
788
789 if( verb ) INFO_message("begin BrainSync [%s]",TIMER) ;
790
791 /* build mask */
792
793 nvox = DSET_NVOX(dsetB) ;
794 mm = DSET_NVALS(dsetB) ;
795
796 if( maskset != NULL ){ /*** = user-supplied mask ***/
797 vmask = THD_makemask( maskset , 0 , 1.0f,0.0f ) ;
798 DSET_unload(maskset) ;
799 if( vmask == NULL )
800 ERROR_exit("Can't make mask from -mask dataset '%s'",DSET_BRIKNAME(maskset)) ;
801 nvmask = THD_countmask( nvox , vmask ) ;
802 if( verb ) ININFO_message("%d voxels in the spatial mask",nvmask) ;
803 if( nvmask == 0 )
804 ERROR_exit("Mask from -mask dataset %s has 0 voxels",DSET_BRIKNAME(maskset)) ;
805 } else { /*** mask = all voxels */
806 vmask = (byte *)malloc(sizeof(byte)*(nvox+2)) ; nvmask = nvox ;
807 for( jj=0 ; jj < nvox ; jj++ ) vmask[jj] = 1 ;
808 if( verb )
809 ININFO_message("no -mask option ==> processing all %d voxels in dataset",nvox);
810 }
811 MEMORY_CHECK("P") ;
812
813 /* find/eliminate (remove from vmask) any voxels with all-constant vectors */
814
815 if( verb > 1 )
816 ININFO_message("looking for all-constant time series [%s]",TIMER) ;
817 bar = (float *)malloc(sizeof(float)*mm) ;
818 for( kk=0 ; kk < nvox ; kk++ ){
819 if( vmask[kk] ){
820 THD_extract_array( kk , dsetB , 0 , bar ) ;
821 if( is_vector_constant(mm,bar) ){ vmask[kk] = 0 ; ncut++ ; continue ; }
822 THD_extract_array( kk , dsetC , 0 , bar ) ;
823 if( is_vector_constant(mm,bar) ){ vmask[kk] = 0 ; ncut++ ; continue ; }
824 }
825 }
826 free(bar) ;
827 if( ncut > 0 ){
828 nvmask = THD_countmask(nvox,vmask) ;
829 if( verb )
830 ININFO_message("removed %d voxel%s for being constant in time" ,
831 ncut , (ncut > 1) ? "s" : "\0" ) ;
832 }
833 if( nvmask <= 2*mm )
834 ERROR_exit("%d is not enough voxels in mask to process :(",nvmask) ;
835
836 /* load datasets into matrices */
837
838 if( verb > 1 )
839 ININFO_message("loading datasets into matrices [%s]",TIMER) ;
840
841 bmat = (float *)calloc( sizeof(float) , mm*nvmask ) ;
842 cmat = (float *)calloc( sizeof(float) , mm*nvox ) ; /* extra big */
843 qmat = (float *)calloc( sizeof(float) , mm*mm ) ;
844
845 for( ii=0 ; ii < mm ; ii++ ){
846 bim = THD_extract_float_brick( ii , dsetB ) ; bar = MRI_FLOAT_PTR(bim) ;
847 cim = THD_extract_float_brick( ii , dsetC ) ; car = MRI_FLOAT_PTR(cim) ;
848 for( kk=jj=0 ; jj < nvox ; jj++ ){
849 if( vmask[jj] ){
850 bmat[ii+kk*mm] = bar[jj] ;
851 cmat[ii+kk*mm] = car[jj] ; kk++ ;
852 }
853 }
854 mri_free(bim) ; mri_free(cim) ;
855 }
856 MEMORY_CHECK("Q") ;
857 DSET_unload(dsetB) ;
858 MEMORY_CHECK("R") ;
859
860 #ifdef DO_DESPIKE /*-----------------------------------------------------*/
861 /* despike matrix columns */
862
863 if( verb > 1 )
864 ININFO_message("despiking time series [%s]",TIMER) ;
865
866 for( ncut=kk=0 ; kk < nvmask ; kk++ ){
867 ncut += despike9( mm , bmat+kk*mm ) ;
868 ncut += despike9( mm , cmat+kk*mm ) ;
869 }
870 if( ncut > 0 )
871 if( verb > 1 )
872 ININFO_message(" squashed %d spike%s from data vectors [%s]",
873 ncut , (ncut==1)?"\0":"s" , TIMER ) ;
874 #endif /* DO_DESPIKE */ /*-----------------------------------------------*/
875
876 /*---------------------------------------------------------------------*/
877 /*-------- compute transform matrix qmat and permutation bperm --------*/
878 /*-------- (that is, do all the work we came here for) --------*/
879 /*---------------------------------------------------------------------*/
880
881 compute_joshi_matrix( mm , nvmask , bmat , cmat , qmat ) ;
882 MEMORY_CHECK("S") ;
883
884 free(bmat) ; /* not needed no more */
885
886 if( verb ) ININFO_message("Brainsync done - computing outputs [%s]",TIMER) ;
887
888 if( Qprefix != NULL ){ /* compute ortset = Qprefix output dataset */
889
890 /* reload input matrix C with ALL voxels, not just mask */
891 if( nvmask < nvox ){
892 if( verb > 1 )
893 ININFO_message("reloading C matrix with all voxels [%s]",TIMER) ;
894 for( ii=0 ; ii < mm ; ii++ ){
895 cim = THD_extract_float_brick( ii , dsetC ) ; car = MRI_FLOAT_PTR(cim) ;
896 for( jj=0 ; jj < nvox ; jj++ ) cmat[ii+jj*mm] = car[jj] ;
897 mri_free(cim) ;
898 }
899 }
900
901 if( verb > 1 )
902 ININFO_message("Q transforming dataset [%s]",TIMER) ;
903
904 cvec = (float *)calloc( sizeof(float) , mm ) ;
905
906 #undef QQ
907 #undef C
908 #define QQ(i,j) qmat[(i)+(j)*mm]
909 #define C(i,j) cmat[(i)+(j)*mm]
910
911 for( jj=0 ; jj < nvox ; jj++ ){
912 for( kk=0 ; kk < mm && C(kk,jj)==0.0f ; kk++ ) ; /*nada*/
913 if( kk == mm ) continue ; /* skip all zero vector */
914 for( ii=0 ; ii < mm ; ii++ ){
915 for( csum=0.0f,kk=0 ; kk < mm ; kk++ ) csum += QQ(ii,kk) * C(kk,jj) ;
916 cvec[ii] = csum ;
917 }
918 if( do_norm ){ /* -normalize column cvec */
919 for( csum=0.0f,ii=0 ; ii < mm ; ii++ ) csum += cvec[ii]*cvec[ii] ;
920 if( csum > 0.0f ){
921 csum = 1.0f/sqrtf(csum) ;
922 for( ii=0 ; ii < mm ; ii++ ) cvec[ii] *= csum ;
923 }
924 }
925 for( ii=0 ; ii < mm ; ii++ ) C(ii,jj) = cvec[ii] ; /* overload matrix */
926 }
927 MEMORY_CHECK("T") ;
928
929 free(cvec) ; free(qmat) ;
930
931 if( verb > 1 )
932 ININFO_message("filling output dataset [%s]",TIMER) ;
933
934 ortset = EDIT_empty_copy( dsetB ) ;
935 EDIT_dset_items( ortset ,
936 ADN_prefix , Qprefix ,
937 ADN_datum_all , MRI_float ,
938 ADN_brick_fac , NULL ,
939 ADN_none ) ;
940 for( ii=0 ; ii < mm ; ii++ ){
941 car = (float *)calloc(sizeof(float),nvox) ;
942 for( jj=0 ; jj < nvox ; jj++ ) car[jj] = C(ii,jj) ;
943 EDIT_substitute_brick( ortset , ii , MRI_float , car ) ;
944 }
945 MEMORY_CHECK("U") ;
946
947 } /*----- end of making ortset -----*/
948
949 /*------------ do it again with the permutation ------------------------*/
950
951 if( Pprefix != NULL && bperm != NULL ){
952 if( verb > 1 )
953 ININFO_message("reloading C matrix with all voxels [%s]",TIMER) ;
954 DSET_load(dsetC) ;
955 for( ii=0 ; ii < mm ; ii++ ){
956 cim = THD_extract_float_brick( ii , dsetC ) ; car = MRI_FLOAT_PTR(cim) ;
957 for( jj=0 ; jj < nvox ; jj++ ) cmat[ii+jj*mm] = car[jj] ;
958 mri_free(cim) ;
959 }
960 cvec = (float *)calloc( sizeof(float) , mm ) ;
961 if( verb > 1 )
962 ININFO_message("permuting C matrix [%s]",TIMER) ;
963 for( jj=0 ; jj < nvox ; jj++ ){
964 for( kk=0 ; kk < mm && C(kk,jj)==0.0f ; kk++ ) ; /*nada*/
965 if( kk == mm ) continue ; /* skip all zero vector */
966 for( ii=0 ; ii < mm ; ii++ ) cvec[ii] = C(bperm[ii],jj) ;
967 if( do_norm ){ /* -normalize */
968 for( csum=0.0f,ii=0 ; ii < mm ; ii++ ) csum += cvec[ii]*cvec[ii] ;
969 if( csum > 0.0f ){
970 csum = 1.0f/sqrtf(csum) ;
971 for( ii=0 ; ii < mm ; ii++ ) cvec[ii] *= csum ;
972 }
973 }
974 for( ii=0 ; ii < mm ; ii++ ) C(ii,jj) = cvec[ii] ;
975 }
976 free(cvec) ; free(bperm) ;
977 if( verb > 1 )
978 ININFO_message("filling permuted dataset [%s]",TIMER) ;
979 permset = EDIT_empty_copy( dsetB ) ;
980 EDIT_dset_items( permset ,
981 ADN_prefix , Pprefix ,
982 ADN_datum_all , MRI_float ,
983 ADN_brick_fac , NULL ,
984 ADN_none ) ;
985 for( ii=0 ; ii < mm ; ii++ ){
986 car = (float *)calloc(sizeof(float),nvox) ;
987 for( jj=0 ; jj < nvox ; jj++ ) car[jj] = C(ii,jj) ;
988 EDIT_substitute_brick( permset , ii , MRI_float , car ) ;
989 }
990 } /*--------------------------------------------------------------*/
991
992 /* datasets are created, so return */
993
994 DSET_unload(dsetC) ; free(cmat) ;
995
996 if( verb ) ININFO_message("outputs computed [%s]",TIMER) ;
997
998 EXRETURN ;
999 }
1000
1001 /*----------------------------------------------------------------------------*/
1002
main(int argc,char * argv[])1003 int main( int argc , char *argv[] )
1004 {
1005 int iarg=1 , mm,nn , nbad=0 ;
1006
1007 /*----------*/
1008
1009 if( argc < 2 || strcasecmp(argv[1],"-help") == 0 )
1010 BSY_help_the_pitiful_user() ;
1011
1012 /*----- bureaucracy -----*/
1013
1014 mainENTRY("3dBrainSync"); machdep();
1015 AFNI_logger("3dBrainSync",argc,argv);
1016 PRINT_VERSION("3dBrainSync"); AUTHOR("Cox the ChronoSynclastic") ;
1017 ct = NI_clock_time() ;
1018 #ifdef USING_MCW_MALLOC
1019 enable_mcw_malloc() ;
1020 #endif
1021
1022 /*----- scan options -----*/
1023
1024 while( iarg < argc && argv[iarg][0] == '-' ){
1025
1026 /*-----*/
1027
1028 if( strcasecmp(argv[iarg],"-quiet") == 0 ){
1029 verb = 0 ; iarg++ ; continue ;
1030 }
1031 if( strcasecmp(argv[iarg],"-verb") == 0 ){
1032 verb++ ; iarg++ ; continue ;
1033 }
1034
1035 /*-----*/
1036
1037 if( strcasecmp(argv[iarg],"-inset1") == 0 || strcasecmp(argv[iarg],"-input1") == 0 ){
1038 if( dsetB != NULL )
1039 ERROR_exit("Can't use option '%s' twice!",argv[iarg]) ;
1040 if( ++iarg >= argc )
1041 ERROR_exit("Need value after option '%s'",argv[iarg-1]) ;
1042
1043 dsetB = THD_open_dataset(argv[iarg]) ;
1044 CHECK_OPEN_ERROR(dsetB,argv[iarg]);
1045 DSET_mallocize(dsetB) ; /* is faster than mmap */
1046 iarg++ ; continue ;
1047 }
1048
1049 /*-----*/
1050
1051 if( strcasecmp(argv[iarg],"-inset2") == 0 || strcasecmp(argv[iarg],"-input2") == 0 ){
1052 if( dsetC != NULL )
1053 ERROR_exit("Can't use option '%s' twice!",argv[iarg]) ;
1054 if( ++iarg >= argc )
1055 ERROR_exit("Need value after option '%s'",argv[iarg-1]) ;
1056
1057 dsetC = THD_open_dataset(argv[iarg]) ;
1058 CHECK_OPEN_ERROR(dsetC,argv[iarg]);
1059 DSET_mallocize(dsetC) ;
1060 iarg++ ; continue ;
1061 }
1062
1063 /*-----*/
1064
1065 if( strcasecmp(argv[iarg],"-mask") == 0 ){
1066 if( maskset != NULL ) ERROR_exit("Can't use option '%s' twice!",argv[iarg]) ;
1067 if( ++iarg >= argc ) ERROR_exit("Need value after option '%s'",argv[iarg-1]) ;
1068 maskset = THD_open_dataset(argv[iarg]) ;
1069 CHECK_OPEN_ERROR(maskset,argv[iarg]) ;
1070 DSET_mallocize(maskset) ;
1071 iarg++ ; continue ;
1072 }
1073
1074 /*-----*/
1075
1076 if( strcasecmp(argv[iarg],"-Qprefix") == 0 ){
1077 char *cpt ;
1078 if( ++iarg >= argc )
1079 ERROR_exit("Need value after option '%s'",argv[iarg-1]) ;
1080 Qprefix = strdup(argv[iarg]) ;
1081 if( !THD_filename_ok(Qprefix) )
1082 ERROR_exit("-Qprefix '%s' is not acceptable :-(",Qprefix) ;
1083 Qbase = strdup(Qprefix) ;
1084 cpt = strstr(Qbase,".nii") ; if( cpt != NULL ) *cpt = '\0' ;
1085 cpt = strstr(Qbase,".1D" ) ; if( cpt != NULL ) *cpt = '\0' ;
1086 iarg++ ; continue ;
1087 }
1088
1089 /*-----*/
1090
1091 if( strcasecmp(argv[iarg],"-Pprefix") == 0 ){
1092 char *cpt ;
1093 if( ++iarg >= argc )
1094 ERROR_exit("Need value after option '%s'",argv[iarg-1]) ;
1095 Pprefix = strdup(argv[iarg]) ;
1096 if( !THD_filename_ok(Pprefix) )
1097 ERROR_exit("-Pprefix '%s' is not acceptable :-(",Pprefix) ;
1098 Pbase = strdup(Pprefix) ;
1099 cpt = strstr(Pbase,".nii") ; if( cpt != NULL ) *cpt = '\0' ;
1100 cpt = strstr(Pbase,".1D" ) ; if( cpt != NULL ) *cpt = '\0' ;
1101 iarg++ ; continue ;
1102 }
1103
1104 /*-----*/
1105
1106 if( strncasecmp(argv[iarg],"-norm",5) == 0 ){
1107 do_norm = 1 ; iarg++ ; continue ;
1108 }
1109
1110 /*--- error! ---*/
1111
1112 ERROR_exit(
1113 "[3dBrainSync] Don't know what to do with option '%s' :(",argv[iarg]) ;
1114
1115 } /* end of option scanning loop */
1116
1117 /*----- error checking the stoopid luser; nbad = fatal error count -----*/
1118
1119 /* no outputs? */
1120
1121 if( Qprefix == NULL && Pprefix == NULL ){
1122 ERROR_message("no -Qprefix or -Pprefix option? Nothing to compute!") ;
1123 nbad++ ;
1124 }
1125
1126 /* bad output names? */
1127
1128 if( Qprefix != NULL && Pprefix != NULL ){
1129 if( strcmp(Qprefix,Pprefix) == 0 ){
1130 ERROR_message("-Qprefix and -Pprefix cannot be the same!") ;
1131 nbad++ ;
1132 } else if( strcmp(Qbase,Pbase) == 0 ){
1133 ERROR_message("-Qprefix and -Pprefix are too similar!") ;
1134 nbad++ ;
1135 }
1136 }
1137
1138 /* no inputs? */
1139
1140 if( dsetB == NULL ){ ERROR_message("no -inset1 dataset?") ; nbad++ ; }
1141 if( dsetC == NULL ){ ERROR_message("no -inset2 dataset?") ; nbad++ ; }
1142
1143 /* grids don't match */
1144
1145 if( !EQUIV_GRIDXYZ(dsetB,dsetC) ){
1146 ERROR_message("-inset1 and -inset2 datasets are not on same 3D grid :(") ;
1147 nbad++ ;
1148 }
1149
1150 /* number of time points is bad? */
1151
1152 mm = DSET_NVALS(dsetB) ;
1153 if( mm < 19 ){
1154 ERROR_message("-inset1 only has %d time points -- need at least 19 :(",mm);
1155 nbad++ ;
1156 }
1157 if( DSET_NVALS(dsetC) != mm ){
1158 ERROR_message(
1159 "-inset1 has %d time points; -inset2 has %d -- they should match :(",
1160 mm , DSET_NVALS(dsetC) ) ;
1161 nbad++ ;
1162 }
1163
1164 /* mask dataset grid doesn't match? */
1165
1166 if( maskset != NULL && !EQUIV_GRIDXYZ(dsetB,maskset) ){
1167 ERROR_message("-mask and -inset1 datsets are NOT on the same 3D grid :(") ;
1168 nbad++ ;
1169 }
1170
1171 /* any fatal errors == DEATH */
1172
1173 if( nbad > 0 ) ERROR_exit("Can't continue after above ERROR%s" ,
1174 (nbad==1) ? "\0" : "s" ) ;
1175
1176 /*----- Load datasets -----*/
1177
1178 if( verb )
1179 INFO_message("start reading datasets [%s]",TIMER) ;
1180 DSET_load(dsetB) ; CHECK_LOAD_ERROR(dsetB) ;
1181 DSET_load(dsetC) ; CHECK_LOAD_ERROR(dsetC) ;
1182 if( maskset != NULL ){ DSET_load(maskset) ; CHECK_LOAD_ERROR(maskset) ; }
1183
1184 /*----- process the data, get output datasets (all the work ) -----*/
1185
1186 BSY_process_data() ;
1187
1188 /* didn't get either output dataset? Weird */
1189
1190 if( ortset == NULL && permset == NULL )
1191 ERROR_exit("Processing the data failed for some reason :(") ;
1192
1193 if( verb > 1 )
1194 ININFO_message("Writing results [%s]" , TIMER ) ;
1195
1196 if( ortset != NULL ){
1197 tross_Copy_History( dsetB , ortset ) ;
1198 tross_Make_History( "3dBrainSync" , argc,argv , ortset ) ;
1199 DSET_write(ortset) ;
1200 if( verb ) WROTE_DSET(ortset) ;
1201 }
1202
1203 if( permset != NULL ){
1204 tross_Copy_History( dsetB , permset ) ;
1205 tross_Make_History( "3dBrainSync" , argc,argv , permset ) ;
1206 DSET_write(permset) ;
1207 if( verb ) WROTE_DSET(permset) ;
1208 }
1209
1210 /* take a victory lap */
1211
1212 if( verb )
1213 INFO_message("=== 3dBrainSync clock time = %s" , TIMER ) ;
1214
1215 exit(0) ;
1216 }
1217