1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4 #include <unistd.h>
5 #include <debugtrace.h>
6 #include <gsl/gsl_randist.h>
7 #include <gsl/gsl_rng.h>
8 #include "mrilib.h" // AFNIadd
9 #include "3ddata.h" // AFNIadd
10 #include "TrackIO.h"
11 #include "DoTrackit.h"
12
13
14
15 /*
16 Added in Feb,2015.
17
18 Trim away record of bundles that are 'too small', whose threshold
19 can be defined by the user.
20
21 */
ByeByeBundle(int A,int B,int NET,int ** Prob_grid,float *** Prob_grid_L,float *** Param_grid,int L_paramgrid,int *** NETROI,int L_netroi,int * NROI)22 int ByeByeBundle( int A,
23 int B,
24 int NET,
25 int **Prob_grid,
26 float ***Prob_grid_L,
27 float ***Param_grid,
28 int L_paramgrid,
29 int ***NETROI,
30 int L_netroi,
31 int *NROI )
32 {
33 int i, lll;
34
35 // magical matrix location
36 lll = MatrInd_to_FlatUHT(A, B, NROI[NET]);
37
38 Prob_grid[NET][lll] = 0;
39
40 Prob_grid_L[NET][lll][0] = 0; // Nov,2016
41 Prob_grid_L[NET][lll][1] = 0;
42
43 for( i=0 ; i<L_paramgrid ; i++ )
44 Param_grid[NET][lll][i] = 0;
45
46 for( i=0 ; i<L_netroi ; i++ )
47 if( NETROI[i][NET][lll] )
48 NETROI[i][NET][lll] = 0;
49
50 RETURN(1);
51 }
52
53
54 /*
55 Order no longer matters here for 'U'HT.
56 We check for the order internally
57 Inputs: NxN square matr coor -> (ROW, COL, N)
58 Output: flat UHT index value PLUS 1
59 *This one goes along the diagonal first!*
60 Useful for unique enumeration of tracts, hence the plus one
61 */
MatrInd_to_FlatUHT_DIAG_P1(int i,int j,int N)62 int MatrInd_to_FlatUHT_DIAG_P1(int i, int j, int N)
63 {
64 int lll;
65
66 if( i<=j ) {
67 i = j-i;
68 lll = i*N+j; // sq matr coor
69 lll-= (i*(i+1))/2; // fix for tridiag.
70 }
71 else {
72 j = i-j;
73 lll = j*N+i; // sq matr coor
74 lll-= (j*(j+1))/2; // fix for tridiag.
75 }
76 lll+=1;
77
78 return lll;
79 }
80
81 /*
82 For when there are multiple connections-- take the unique UHT decomp
83 from P1, and then add N(N-1)/2 to it. Specifically when the max num
84 of overlaps is 2, but we'll also have multi PLUS two of the overlaps
85 in the string label.
86 Max value to be output (even after the +1): N*N.
87 */
MatrInd_to_FlatUHT_DIAG_M(int i,int j,int N)88 int MatrInd_to_FlatUHT_DIAG_M(int i, int j, int N)
89 {
90 int lll;
91
92 lll = MatrInd_to_FlatUHT_DIAG_P1(i,j,N);
93 // add N(N-1)/2 to that value.
94 lll+= FlatUHT_Len(N);
95 lll-= N;
96
97 return lll;
98 }
99
100 /*
101 Order no longer matters here for 'U'HT.
102 We check for the order internally
103 Inputs: NxN square matr coor -> (ROW, COL, N)
104 Output: flat UHT index value
105 */
MatrInd_to_FlatUHT(int i,int j,int N)106 int MatrInd_to_FlatUHT(int i, int j, int N)
107 {
108 int lll;
109
110 if( i<=j ) {
111 lll = i*N+j; // sq matr coor
112 lll-= (i*(i+1))/2; // fix for tridiag.
113 }
114 else {
115 lll = j*N+i; // sq matr coor
116 lll-= (j*(j+1))/2; // fix for tridiag.
117 }
118
119 return lll;
120 }
121
FlatUHT_Len(int N)122 int FlatUHT_Len(int N)
123 {
124 int out;
125 out = N*(N+1); // out must have been an even number...
126 out/= 2;
127 return out;
128 }
129
130
131 // Some funniness to deal with possibility of having biggest label of
132 // a tractogr ROI be >M, where M is actual number of ROIs. Before, we
133 // required labels to be 1..M, but now we'll free that restriction up.
134 // REF is input data set of ROI labels
135 // NUMROI ends with number of ROIs per brik
136 // - dimensions: N_briks x 1
137 // - stores value of M per brik
138 // ROILIST ends with having ordered list of ROILABEL ints
139 // - dimensions: N_briks x (max ROI label)+1
140 // - data dims: N_briks x M+1
141 // (where the +1s are because we go 1....value)
142 // INVLIST has the ith values at the actual input locations... hence gives
143 // inverse info to ROILIST
144 // INVROI starts with max val per brik; ends with keeping
145 // track of max input index per brik
146 //
147 // for now, all negative ROI indexing is handled in using MAPROI
148 // don't think negative integers can be in refset for 3dROIMaker yet, though
ViveLeRoi(THD_3dim_dataset * REF,int ** ROILIST,int ** INVLIST,int * NUMROI,int * INVROI)149 int ViveLeRoi(THD_3dim_dataset *REF, int **ROILIST, int **INVLIST,
150 int *NUMROI, int *INVROI)
151 { // NB what index boundaries are for various loops here...
152 int Nbrik=0,Nvox=0;
153 int i,j,k,m;
154 int nrois=0;
155
156 Nbrik = DSET_NVALS(REF);
157 Nvox = DSET_NVOX(REF);
158
159
160 for( m=0 ; m<Nbrik ; m++ )
161 for( i=0 ; i<Nvox ; i++ )
162 if( THD_get_voxel(REF,i,m) > 0.5 ) {
163 ROILIST[m][(int) THD_get_voxel(REF,i,m)]=1;
164 }
165
166 // So, all M rois per brik are now marked. Go through them,
167 // condense them into a list by array indices 0..M-1. NB, in this
168 // case, the actual *value* of the ROI label might be >M.
169 // The value of M is then stored for real in NUMROI.
170 for( m=0 ; m<Nbrik ; m++ ) {
171 j=1;
172 for( i=1 ; i<=INVROI[m] ; i++ ){
173 if( ROILIST[m][i] == 1 ){
174 ROILIST[m][j]=i;
175 INVLIST[m][i]=j;
176 j++;
177 }
178 }
179
180 if(INVROI[m] < j-1) // should never happen...
181 ERROR_exit("Problem with ROI labels! Badness in reading/counting.");
182
183 // reset value with real total number, not just max label
184 NUMROI[m] = j-1;
185 }
186
187 RETURN(1);
188 };
189
190
191 /*
192 check a given voxel about whether the track running through it is in
193 an entered NOT-mask
194 */
CheckNotMask(int id,int br,short ** amask,int AO)195 int CheckNotMask(int id, int br, short **amask, int AO){
196 int out=0;
197
198 if(AO){ // N_nets==NOTMASK
199 if(amask[id][br])
200 out=1;
201 }
202 else // single NOTMASK for all
203 if(amask[id][0])
204 out=1;
205
206 return out;
207 }
208
209
210 /*
211 varied length version
212
213 */
ScoreTrackGrid_M(float *** PG,int idx,int h,int C,THD_3dim_dataset ** inset,int bot,int top)214 int ScoreTrackGrid_M( float ***PG, int idx, int h, int C,
215 THD_3dim_dataset **inset, int bot, int top)
216 {
217 int i,j=1;
218 float READS_fl=0;
219
220 PG[h][C][0]+= 1.0; // N_vox
221
222 for( i=bot ; i<top ; i++ ) {
223
224 //mu and std of FA,MD,RD,L1
225 PG[h][C][j]+=
226 THD_get_voxel(inset[i],idx,0);
227 PG[h][C][j+1]+=
228 (float) pow(THD_get_voxel(inset[i],idx,0),2);
229 j+=2;
230 }
231
232 RETURN(0);
233 }
234
235
236
237 /* OLD
238 int ScoreTrackGrid_M( float ****PG,int idx, int h, int C, int B,
239 THD_3dim_dataset **inset, int bot, int top,
240 int DTI)
241 {
242 int i,j=1;
243 float READS_fl=0;
244
245 PG[h][C][B][0]+= 1.0; // N_vox
246
247 for( i=bot ; i<top ; i++ ) {
248
249 //mu and std of FA,MD,RD,L1
250 PG[h][C][B][j]+=
251 THD_get_voxel(inset[i],idx,0);
252 PG[h][C][B][j+1]+=
253 (float) pow(THD_get_voxel(inset[i],idx,0),2);
254 j+=2;
255 }
256
257 if(DTI) { // RD = 0.5 * ( 3*MD - L1 )
258 READS_fl = 0.5*(3.0*THD_get_voxel(inset[i-2],idx,0)-
259 THD_get_voxel(inset[i-1],idx,0));
260 PG[h][C][B][j]+=
261 READS_fl;
262 PG[h][C][B][j+1]+=
263 (float) pow(READS_fl,2);
264 }
265
266 RETURN(1);
267 }*/
268
269
270
271 // *******************************************************************
272
273 // basic format for writing out tracking results of 3dTrackID:
274 // + a file of individual ROI intercepts
275 // + a file of the pairwise combinations
WriteBasicProbFiles(int N_nets,int Ndata,int Nvox,char * prefix,THD_3dim_dataset * insetFA,int * TV_switch,char * voxel_order,int * NROI,int *** NETROI,int *** mskd,int *** INDEX2,int * Dim,THD_3dim_dataset * dsetn,int argc,char * argv[],char *** ROI_STR_LABS,int NameLabelsOut,Dtable * roi_table,int ** roi_labs,int PAIR_POWERON,int NIFTI_OUT)276 int WriteBasicProbFiles(int N_nets, int Ndata, int Nvox,
277 char *prefix,THD_3dim_dataset *insetFA,
278 int *TV_switch,char *voxel_order,int *NROI,
279 int ***NETROI,int ***mskd,int ***INDEX2,int *Dim,
280 THD_3dim_dataset *dsetn,int argc, char *argv[],
281 char ***ROI_STR_LABS, int NameLabelsOut,
282 Dtable *roi_table,
283 int **roi_labs, int PAIR_POWERON,
284 int NIFTI_OUT)
285 {
286
287 int i,j,k,bb,hh,kk,rr,idx;
288 char **prefix_netmap=NULL;
289 char **prefix_netmap2=NULL;
290 char **prefix_dtable=NULL;
291 char *Dtable_str=NULL;
292 Dtable *new_dt=NULL;
293 char mini[50];
294 THD_3dim_dataset *networkMAPS=NULL,*networkMAPS2=NULL;
295 char bric_labs[THD_MAX_NAME];
296 int idx3;
297 FILE *fout1;
298 int MAXOVERLAP = 2;
299 int i1,i2;
300 char EleNameStr[128];
301 int maxROInum, MAXOVERLAP_VAL;
302 int MULTI_ROI = 0;
303
304 // ****** alloc'ing
305 prefix_netmap = calloc( N_nets, sizeof(prefix_netmap));
306 for(i=0 ; i<N_nets ; i++)
307 prefix_netmap[i] = calloc( THD_MAX_NAME, sizeof(char));
308 prefix_netmap2 = calloc( N_nets, sizeof(prefix_netmap2));
309 for(i=0 ; i<N_nets ; i++)
310 prefix_netmap2[i] = calloc( THD_MAX_NAME, sizeof(char));
311 prefix_dtable = calloc( N_nets, sizeof(prefix_dtable));
312 for(i=0 ; i<N_nets ; i++)
313 prefix_dtable[i] = calloc( THD_MAX_NAME, sizeof(char));
314
315 if( ( prefix_netmap== NULL) || ( prefix_netmap2== NULL)
316 || ( prefix_dtable== NULL) ) {
317 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
318 exit(122);
319 }
320
321 // ****** calc/do
322 for( hh=0 ; hh<N_nets ; hh++) {
323
324 // Oct. 2014: if only 1 ROI in set, then just output 0th brick
325 // (because others add no information)
326 // MULTI_ROI acts as both a switch and a number
327 MULTI_ROI = ( NROI[hh] > 1 ) ? 1 : 0;
328
329 if( NIFTI_OUT )
330 sprintf(prefix_netmap[hh],"%s_%03d_PAIRMAP.nii.gz",prefix,hh);
331 else
332 sprintf(prefix_netmap[hh],"%s_%03d_PAIRMAP",prefix,hh);
333
334 // Sept 2014
335 sprintf(prefix_dtable[hh],"%s_%03d_PAIRMAP.niml.lt",prefix, hh);
336 if( roi_table ) {
337 // copy dtable
338 Dtable_str = Dtable_to_nimlstring(roi_table, "VALUE_LABEL_DTABLE");
339 new_dt = Dtable_from_nimlstring(Dtable_str);
340 free(Dtable_str); Dtable_str=NULL;
341 }
342 else {
343 new_dt = new_Dtable( NROI[hh] );
344 }
345 // from either above cases, we have either a non-empty or empty
346 // starter; in either case, check through all our ROIs to add in
347 // more as necessary.
348 for( bb=1 ; bb<=NROI[hh] ; bb++) {
349 snprintf(mini, 50, "%d", roi_labs[hh][bb]);
350 if(!(findin_Dtable_a( mini, new_dt ))) {
351 addto_Dtable(mini, ROI_STR_LABS[hh][bb], new_dt );
352 }
353 }
354
355
356 // just get one of right dimensions!
357 networkMAPS = EDIT_empty_copy( insetFA ) ;
358 // Oct. 2014: if only 1 ROI in set, then just output 0th brick
359 // (because others add no information)
360 if( MULTI_ROI )
361 EDIT_add_bricklist(networkMAPS ,
362 NROI[hh], NULL , NULL , NULL );
363 if( PAIR_POWERON )
364 EDIT_dset_items(networkMAPS,
365 ADN_datum_all , MRI_float ,
366 ADN_none ) ;
367 else
368 EDIT_dset_items(networkMAPS,
369 ADN_datum_all , MRI_short ,
370 ADN_none ) ;
371
372 if( NIFTI_OUT )
373 sprintf(prefix_netmap2[hh],"%s_%03d_INDIMAP.nii.gz",prefix,hh);
374 else
375 sprintf(prefix_netmap2[hh],"%s_%03d_INDIMAP",prefix,hh);
376
377 // just get one of right dimensions!
378 networkMAPS2 = EDIT_empty_copy( insetFA ) ;
379 // Oct. 2014: if only 1 ROI in set, then just output 0th brick
380 // (because others add no information)
381 if( MULTI_ROI )
382 EDIT_add_bricklist(networkMAPS2 ,
383 NROI[hh], NULL , NULL , NULL );
384 EDIT_dset_items(networkMAPS2,
385 ADN_datum_all , MRI_float ,
386 ADN_none ) ;
387
388 float **temp_arrFL=NULL;
389 short int **temp_arrSH=NULL;
390 float **temp_arr2=NULL;
391
392 int **intersec=NULL;
393
394 // first array for all tracks, 2nd for paired ones.
395 // still just need one set of matrices output
396 if( MULTI_ROI ) {
397 if( PAIR_POWERON ) {
398 temp_arrFL = calloc( (NROI[hh]+MULTI_ROI), sizeof(temp_arrFL));
399 for(i=0 ; i<(NROI[hh]+MULTI_ROI) ; i++)
400 temp_arrFL[i] = calloc( Nvox, sizeof(float));
401 }
402 else {
403 temp_arrSH = calloc( (NROI[hh]+MULTI_ROI), sizeof(temp_arrSH));
404 for(i=0 ; i<(NROI[hh]+MULTI_ROI) ; i++)
405 temp_arrSH[i] = calloc( Nvox, sizeof(short int));
406 }
407 }
408
409 temp_arr2 = calloc( (NROI[hh]+MULTI_ROI), sizeof(temp_arr2));
410 for(i=0 ; i<(NROI[hh]+MULTI_ROI) ; i++)
411 temp_arr2[i] = calloc( Nvox, sizeof(float));
412
413 if( temp_arr2 == NULL ) {
414 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
415 exit(122);
416 }
417
418 // use this per vox
419 if( MULTI_ROI ){
420 intersec = calloc( (NROI[hh]+MULTI_ROI), sizeof(intersec));
421 for(i=0 ; i<(NROI[hh]+MULTI_ROI) ; i++)
422 intersec[i] = calloc(MAXOVERLAP+1, sizeof( int ));
423
424 if( PAIR_POWERON ) {
425 if( temp_arrFL == NULL ) {
426 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
427 exit(122);
428 }
429 }
430 else {
431 if( temp_arrSH == NULL ) {
432 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
433 exit(122);
434 }
435 }
436 if( intersec == NULL ) {
437 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
438 exit(122);
439 }
440 }
441
442
443 // for string label outputs
444 maxROInum = roi_labs[hh][NROI[hh]];
445 // should never happen normally -- even with _M_<->A<->B
446 MAXOVERLAP_VAL = maxROInum*maxROInum + 1;
447
448 idx=0;
449 for( k=0 ; k<Dim[2] ; k++ )
450 for( j=0 ; j<Dim[1] ; j++ )
451 for( i=0 ; i<Dim[0] ; i++ ) {
452 // allow for more than one `connector' tract
453 if(mskd[i][j][k]) {
454 for( bb=1 ; bb<=NROI[hh] ; bb++)
455 for( rr=bb ; rr<=NROI[hh] ; rr++) {
456 idx3 = MatrInd_to_FlatUHT(bb-1,rr-1,NROI[hh]);
457 if(NETROI[INDEX2[i][j][k]][hh][idx3]>0) {
458 // store connectors
459 if( MULTI_ROI ){
460 if(bb != rr){
461 if(PAIR_POWERON)
462 temp_arrFL[0][idx] = 1.;
463 else
464 temp_arrSH[0][idx] = (short) 1;
465 }
466 }
467
468 // store tracks through any ROI
469 temp_arr2[0][idx] = (float)
470 NETROI[INDEX2[i][j][k]][hh][idx3];
471
472
473 if( MULTI_ROI ) {
474 // then add value if overlap
475 if(bb != rr) {
476 if( PAIR_POWERON ) {// OLD
477 // unique
478 temp_arrFL[bb][idx]+=(float) pow(2,rr);
479 temp_arrFL[rr][idx]+=(float) pow(2,bb);
480 }
481 else{
482 temp_arrSH[bb][idx]+= (short) roi_labs[hh][rr];
483 temp_arrSH[rr][idx]+= (short) roi_labs[hh][bb];
484 intersec[bb][0]+= 1;
485 intersec[rr][0]+= 1;
486 // keep track of which ones get hit,
487 // less than the value
488 if( intersec[bb][0]<=MAXOVERLAP )
489 intersec[bb][intersec[bb][0]] = rr;
490 if( intersec[rr][0]<=MAXOVERLAP )
491 intersec[rr][intersec[rr][0]] = bb;
492 }
493 }
494
495
496 if(bb != rr) {
497 temp_arr2[bb][idx] = (float)
498 NETROI[INDEX2[i][j][k]][hh][idx3];
499 temp_arr2[rr][idx] = (float)
500 NETROI[INDEX2[i][j][k]][hh][idx3];
501 }
502 else
503 temp_arr2[bb][idx] = (float)
504 NETROI[INDEX2[i][j][k]][hh][idx3];
505 }
506
507 }
508 }
509 // continuation of labelling functions for this voxel
510 if( MULTI_ROI )
511 for( bb=1 ; bb<=NROI[hh] ; bb++) {
512
513 if( intersec[bb][0] == 1 ) { // already tabled, ->erase
514 intersec[bb][0] = intersec[bb][1] = 0;
515 }
516 else if ( intersec[bb][0] == 2 ){ // check labeltable
517
518 i1 = roi_labs[hh][intersec[bb][1]];
519 i2 = roi_labs[hh][intersec[bb][2]];
520 temp_arrSH[bb][idx] =
521 MatrInd_to_FlatUHT_DIAG_P1( i1-1,
522 i2-1,
523 maxROInum);
524 snprintf(mini, 50, "%d", (int) temp_arrSH[bb][idx]);
525
526 if(!(findin_Dtable_a( mini, new_dt ))) {
527 snprintf( EleNameStr, 128, "%s<->%s",
528 ROI_STR_LABS[hh][intersec[bb][1]],
529 ROI_STR_LABS[hh][intersec[bb][2]]);
530 addto_Dtable(mini, EleNameStr, new_dt );
531 }
532 intersec[bb][0]=intersec[bb][1]=intersec[bb][2]=0;
533
534 }
535 else if( intersec[bb][0] ) { // 'multi' + 2 label names
536
537 i1 = roi_labs[hh][intersec[bb][1]];
538 i2 = roi_labs[hh][intersec[bb][2]];
539 // use the multi function now.
540 temp_arrSH[bb][idx] =
541 MatrInd_to_FlatUHT_DIAG_M( i1-1,
542 i2-1,
543 maxROInum);
544 snprintf(mini, 50, "%d", (int) temp_arrSH[bb][idx]);
545
546 if(!(findin_Dtable_a( mini, new_dt ))) {
547 snprintf( EleNameStr, 128, "%s<->%s<->%s",
548 "_M_",
549 ROI_STR_LABS[hh][intersec[bb][1]],
550 ROI_STR_LABS[hh][intersec[bb][2]]);
551 addto_Dtable(mini, EleNameStr, new_dt );
552 }
553 intersec[bb][0]=intersec[bb][1]=intersec[bb][2]=0;
554
555 }
556
557 }
558 }
559 idx+=1;
560 }
561
562
563 // FIRST THE PAIR CONNECTORS
564 if ( MULTI_ROI ) {
565 if( PAIR_POWERON ) {// OLD
566 EDIT_substitute_brick(networkMAPS, 0, MRI_float, temp_arrFL[0]);
567 temp_arrFL[0]=NULL;
568 }
569 else {
570 EDIT_substitute_brick(networkMAPS, 0, MRI_short, temp_arrSH[0]);
571 temp_arrSH[0]=NULL;
572 }
573 }
574
575 // THEN THE INDIVID TRACKS
576 EDIT_substitute_brick(networkMAPS2, 0, MRI_float, temp_arr2[0]);
577 temp_arr2[0]=NULL;
578
579 if( MULTI_ROI )
580 for( bb=1 ; bb<=NROI[hh] ; bb++) {
581 if( PAIR_POWERON ) {// OLD
582 EDIT_substitute_brick(networkMAPS,bb,MRI_float,temp_arrFL[bb]);
583 temp_arrFL[bb]=NULL; // to not get into trouble...
584 }
585 else {
586 EDIT_substitute_brick(networkMAPS,bb,MRI_short,temp_arrSH[bb]);
587 temp_arrSH[bb]=NULL; // to not get into trouble...
588 }
589
590 // Sept 2014: updating labelling
591 sprintf(bric_labs,"AND_%s",ROI_STR_LABS[hh][bb]);
592 EDIT_BRICK_LABEL(networkMAPS, bb, bric_labs); // labels, PAIR
593
594 EDIT_substitute_brick(networkMAPS2, bb, MRI_float, temp_arr2[bb]);
595 sprintf(bric_labs,"OR_%s",ROI_STR_LABS[hh][bb]);
596 EDIT_BRICK_LABEL(networkMAPS2, bb, bric_labs); // labels, INDI
597
598 temp_arr2[bb]=NULL; // to not get into trouble...
599 }
600
601 if(TV_switch[0] || TV_switch[1] || TV_switch[2]) {
602 if( MULTI_ROI ){
603 dsetn = r_new_resam_dset(networkMAPS, NULL, 0.0, 0.0, 0.0,
604 voxel_order, RESAM_NN_TYPE,
605 NULL, 1, 0);
606 DSET_delete(networkMAPS);
607 networkMAPS=dsetn;
608 dsetn=NULL;
609 }
610
611 dsetn = r_new_resam_dset(networkMAPS2, NULL, 0.0, 0.0, 0.0,
612 voxel_order, RESAM_NN_TYPE,
613 NULL, 1, 0);
614 DSET_delete(networkMAPS2);
615 networkMAPS2=dsetn;
616 dsetn=NULL;
617 }
618
619 if( MULTI_ROI ) {
620 EDIT_dset_items(networkMAPS,
621 ADN_prefix , prefix_netmap[hh] ,
622 ADN_brick_label_one , "AND_all",
623 ADN_none ) ;
624
625 // Sept 2014
626 if( !PAIR_POWERON ) {
627 Dtable_str = Dtable_to_nimlstring(new_dt, "VALUE_LABEL_DTABLE");
628 destroy_Dtable(new_dt); new_dt = NULL;
629 THD_set_string_atr( networkMAPS->dblk ,
630 "VALUE_LABEL_DTABLE" , Dtable_str);
631
632 if( (fout1 = fopen(prefix_dtable[hh], "w")) == NULL) {
633 fprintf(stderr, "Error opening file %s.",prefix_dtable[hh]);
634 exit(19);
635 }
636 fprintf(fout1,"%s",Dtable_str);
637 fclose(fout1);
638
639 free(Dtable_str); Dtable_str = NULL;
640 }
641
642 THD_load_statistics(networkMAPS);
643 if( !THD_ok_overwrite() && THD_is_ondisk(DSET_HEADNAME(networkMAPS)) )
644 ERROR_exit("Can't overwrite existing dataset '%s'",
645 DSET_HEADNAME(networkMAPS));
646 tross_Make_History("3dTrackID", argc, argv, networkMAPS);
647 THD_write_3dim_dataset(NULL, NULL, networkMAPS, True);
648 DSET_delete(networkMAPS);
649
650 if(PAIR_POWERON){
651 for( i=0 ; i<NROI[hh]+MULTI_ROI ; i++) // free all
652 free(temp_arrFL[i]);
653 free(temp_arrFL);
654 }
655 else{
656 for( i=0 ; i<NROI[hh]+MULTI_ROI ; i++) // free all
657 free(temp_arrSH[i]);
658 free(temp_arrSH);
659 }
660 for( i=0 ; i<NROI[hh]+MULTI_ROI ; i++) // free all
661 free(intersec[i]);
662 free(intersec);
663 intersec = NULL;
664 }
665
666 EDIT_dset_items(networkMAPS2,
667 ADN_prefix , prefix_netmap2[hh] ,
668 ADN_brick_label_one , "OR_all",
669 ADN_none ) ;
670 THD_load_statistics(networkMAPS2);
671 if( !THD_ok_overwrite() && THD_is_ondisk(DSET_HEADNAME(networkMAPS2)) )
672 ERROR_exit("Can't overwrite existing dataset '%s'",
673 DSET_HEADNAME(networkMAPS2));
674 tross_Make_History("3dTrackID", argc, argv, networkMAPS2);
675 THD_write_3dim_dataset(NULL, NULL, networkMAPS2, True);
676 DSET_delete(networkMAPS2);
677 for( i=0 ; i<NROI[hh]+MULTI_ROI ; i++) // free all
678 free(temp_arr2[i]);
679 free(temp_arr2);
680
681 }
682
683 // ****** freeing **********
684
685 if(MULTI_ROI){
686 for( i=0 ; i<N_nets ; i++) {
687 free(prefix_netmap[i]);
688 }
689 free(prefix_netmap);
690 free(networkMAPS);
691 }
692 for( i=0 ; i<N_nets ; i++) {
693 free(prefix_dtable[i]); // free in all cases because where initialized
694 free(prefix_netmap2[i]);
695 }
696 free(prefix_netmap2);
697 free(prefix_dtable);
698 free(networkMAPS2);
699
700 RETURN(1);
701 }
702
703 // second format for writing out tracking results of 3dTrackID:
704 // each pairwise map in own file
705 // (will be in new directory with prefix name)
WriteIndivProbFiles(int N_nets,int Ndata,int Nvox,int ** Prob_grid,char * prefix,THD_3dim_dataset * insetFA,int * TV_switch,char * voxel_order,int * NROI,int *** NETROI,int *** mskd,int *** INDEX2,int * Dim,THD_3dim_dataset * dsetn,int argc,char * argv[],float *** Param_grid,int DUMP_TYPE,int DUMP_ORIG_LABS,int ** ROI_LABELS,int POST_IT,char *** ROI_STR_LAB,int NameLabelsOut,int NIFTI_OUT)706 int WriteIndivProbFiles(int N_nets, int Ndata, int Nvox, int **Prob_grid,
707 char *prefix,THD_3dim_dataset *insetFA,
708 int *TV_switch,char *voxel_order,int *NROI,
709 int ***NETROI,int ***mskd,int ***INDEX2,int *Dim,
710 THD_3dim_dataset *dsetn,int argc, char *argv[],
711 float ***Param_grid, int DUMP_TYPE,
712 int DUMP_ORIG_LABS, int **ROI_LABELS, int POST_IT,
713 char ***ROI_STR_LAB, int NameLabelsOut,
714 int NIFTI_OUT)
715 {
716
717 int i,j,k,bb,hh,rr,ii,jj,kk;
718 int idx,idx2,count;
719 char ***prefix_netmap=NULL;
720 char ***prefix_dump=NULL;
721 THD_3dim_dataset *networkMAPS=NULL;
722 int *N_totpair=NULL;
723 int sum_pairs=0;
724 FILE *fout;
725 int idx3;
726 int OUT_FLOAT_MAP = 0;
727 //char str1[128]={""}, str2[128]={""};
728 //char *str_lab1=NULL, *str_lab2=NULL;
729
730
731 if( POST_IT || ( DUMP_TYPE==4) )
732 OUT_FLOAT_MAP = 1;
733
734
735 N_totpair = (int *)calloc(N_nets, sizeof(int));
736 // str_lab1 = (char *)calloc(100, sizeof(char));
737 //str_lab2 = (char *)calloc(100, sizeof(char));
738
739 // find out how many networks we'll be outputting.
740 // this means going through UHT part of probgrid, looking for nonzeros
741 for( k=0 ; k<N_nets ; k++)
742 for( i=0 ; i<NROI[k] ; i++ )
743 for( j=i ; j<NROI[k] ; j++ ) // include diags
744 if(Prob_grid[k][MatrInd_to_FlatUHT(i,j,NROI[k])]>0){
745 N_totpair[k]+=1;
746 sum_pairs+=1;
747 }
748
749 if( sum_pairs==0 )
750 INFO_message("No pairs in any network to output.");
751 else {
752
753 INFO_message("Writing out individual files to: %s/NET_ ...",prefix);
754
755
756 // ****** alloc'ing
757 prefix_netmap = (char ***) calloc( N_nets, sizeof(char **) );
758 for ( i = 0 ; i < N_nets ; i++ )
759 prefix_netmap[i] = (char **) calloc( N_totpair[i], sizeof(char *) );
760 for ( i = 0 ; i < N_nets ; i++ )
761 for ( j = 0 ; j < N_totpair[i] ; j++ )
762 prefix_netmap[i][j] = (char *) calloc( THD_MAX_NAME, sizeof(char) );
763 prefix_dump = (char ***) calloc( N_nets, sizeof(char **) );
764 for ( i = 0 ; i < N_nets ; i++ )
765 prefix_dump[i] = (char **) calloc( N_totpair[i], sizeof(char *) );
766 for ( i = 0 ; i < N_nets ; i++ )
767 for ( j = 0 ; j < N_totpair[i] ; j++ )
768 prefix_dump[i][j] = (char *) calloc( THD_MAX_NAME, sizeof(char) );
769
770 if( ( prefix_netmap== NULL) || ( prefix_dump== NULL)) {
771 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
772 exit(122);
773 }
774
775 mkdir(prefix, 0777);
776 // ****** calc/do, loop through networks
777 for( hh=0 ; hh<N_nets ; hh++) {
778 count=0;
779 for( i=0 ; i<NROI[hh] ; i++ )
780 for( j=i ; j<NROI[hh] ; j++ ) {// include diags
781 idx3 = MatrInd_to_FlatUHT(i,j,NROI[hh]);
782 if(Prob_grid[hh][idx3]>0) {
783 if( ROI_STR_LAB && NameLabelsOut ) {
784 if( NIFTI_OUT )
785 snprintf(prefix_netmap[hh][count], THD_MAX_NAME,
786 "%s/NET_%03d_ROI_%s__%s.nii.gz", prefix, hh,
787 ROI_STR_LAB[hh][i+1], ROI_STR_LAB[hh][j+1]);
788 else
789 snprintf(prefix_netmap[hh][count], THD_MAX_NAME,
790 "%s/NET_%03d_ROI_%s__%s", prefix, hh,
791 ROI_STR_LAB[hh][i+1], ROI_STR_LAB[hh][j+1]);
792 }
793 else if(!DUMP_ORIG_LABS) {
794 if( NIFTI_OUT )
795 snprintf(prefix_netmap[hh][count], THD_MAX_NAME,
796 "%s/NET_%03d_ROI_%03d__%03d.nii.gz",
797 prefix,hh,i+1,j+1);
798 else
799 snprintf(prefix_netmap[hh][count], THD_MAX_NAME,
800 "%s/NET_%03d_ROI_%03d__%03d",
801 prefix,hh,i+1,j+1);
802 }
803 else{
804 if( NIFTI_OUT )
805 snprintf(prefix_netmap[hh][count], THD_MAX_NAME,
806 "%s/NET_%03d_ROI_%03d__%03d.nii.gz",prefix,hh,
807 ROI_LABELS[hh][i+1],ROI_LABELS[hh][j+1]);
808 else
809 snprintf(prefix_netmap[hh][count], THD_MAX_NAME,
810 "%s/NET_%03d_ROI_%03d__%03d",prefix,hh,
811 ROI_LABELS[hh][i+1],ROI_LABELS[hh][j+1]);
812 }
813
814
815
816
817 // single brik, byte map
818 networkMAPS = EDIT_empty_copy( insetFA ) ;
819 EDIT_dset_items(networkMAPS,
820 ADN_datum_all , MRI_byte ,
821 ADN_none ) ;
822
823 sprintf(prefix_dump[hh][count],
824 "%s/NET_%03d_ROI_%03d__%03d.dump",
825 prefix,hh,i+1,j+1);
826
827 float *temp_arr_FL=NULL;
828 byte *temp_arr_BY=NULL;
829
830 // first array for all tracks, 2nd for paired ones.
831 // still just need one set of matrices output
832 if(OUT_FLOAT_MAP) {
833 // will be single brik output
834 temp_arr_FL = (float *)calloc(Nvox, sizeof(float));
835 if( temp_arr_FL== NULL ) {
836 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
837 exit(120);
838 }
839 }
840 else {
841 // will be single brik output
842 temp_arr_BY = (byte *)calloc(Nvox, sizeof(byte));
843 if( temp_arr_BY== NULL ) {
844 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
845 exit(121);
846 }
847 }
848
849 int **temp_arr2=NULL;
850 // we know how many vox per WM ROI, alloc that much
851 // for the to-be-dumped mask
852 temp_arr2=calloc(Param_grid[hh][idx3][0],
853 sizeof(temp_arr2));
854 for(bb=0 ; bb<Param_grid[hh][idx3][0] ; bb++)
855 temp_arr2[bb] = calloc(4, sizeof(int)); //x,y,z,1
856
857 if( temp_arr2 == NULL ) {
858 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
859 exit(122);
860 }
861
862 idx=0;
863 idx2=0;
864 for( kk=0 ; kk<Dim[2] ; kk++ )
865 for( jj=0 ; jj<Dim[1] ; jj++ )
866 for( ii=0 ; ii<Dim[0] ; ii++ ) {
867 if(mskd[ii][jj][kk]) {
868 idx3 = MatrInd_to_FlatUHT(i,j,NROI[hh]);
869 if(NETROI[INDEX2[ii][jj][kk]][hh][idx3]>0) {
870 // store locations
871 if(OUT_FLOAT_MAP)
872 temp_arr_FL[idx] = (float)
873 NETROI[INDEX2[ii][jj][kk]][hh][idx3];
874 else
875 temp_arr_BY[idx] = 1;
876 temp_arr2[idx2][0] = ii;
877 temp_arr2[idx2][1] = jj;
878 temp_arr2[idx2][2] = kk;
879 temp_arr2[idx2][3] = 1;
880 idx2++;
881 }
882 }
883 idx++;
884 }
885
886 if(idx2 != Param_grid[hh][idx3][0])
887 printf("ERROR IN COUNTING! Netw,ROI,ROI= (%d, %d, %d); "
888 "idx2 %d != %d paramgrid.\n",
889 hh,i,j,idx2,(int) Param_grid[hh][idx3][0]);
890
891 if(OUT_FLOAT_MAP){
892 EDIT_substitute_brick(networkMAPS, 0, MRI_float,
893 temp_arr_FL);
894 temp_arr_FL=NULL; // to not get into trouble...
895 }
896 else{
897 EDIT_substitute_brick(networkMAPS, 0, MRI_byte,
898 temp_arr_BY);
899 temp_arr_BY=NULL; // to not get into trouble...
900 }
901
902 if(TV_switch[0] || TV_switch[1] || TV_switch[2]) {
903 dsetn = r_new_resam_dset(networkMAPS, NULL, 0.0, 0.0, 0.0,
904 voxel_order, RESAM_NN_TYPE,
905 NULL, 1, 0);
906 DSET_delete(networkMAPS);
907 networkMAPS=dsetn;
908 dsetn=NULL;
909
910 for( bb=0 ; bb<3 ; bb++ )
911 if(TV_switch[bb])
912 for( rr=0 ; rr<idx2 ; rr++)
913 temp_arr2[rr][bb] = Dim[bb]-1-temp_arr2[rr][bb];
914
915 }
916 EDIT_dset_items(networkMAPS,
917 ADN_prefix , prefix_netmap[hh][count] ,
918 ADN_brick_label_one , "ROI_pair",
919 ADN_none ) ;
920
921 // output AFNI if D_T=2 or 3 -> or now 4
922 if(DUMP_TYPE>1) {
923 THD_load_statistics(networkMAPS);
924 if( !THD_ok_overwrite() &&
925 THD_is_ondisk(DSET_HEADNAME(networkMAPS)) )
926 ERROR_exit("Can't overwrite existing dataset '%s'",
927 DSET_HEADNAME(networkMAPS));
928 tross_Make_History("3dTrackID", argc, argv, networkMAPS);
929 THD_write_3dim_dataset(NULL, NULL, networkMAPS, True);
930 }
931
932 DSET_delete(networkMAPS);
933 if(OUT_FLOAT_MAP)
934 free(temp_arr_FL);
935 else
936 free(temp_arr_BY);
937
938 // THEN THE DUMP TEXT FILE
939 // if D_T = 1 or 3
940 if( (DUMP_TYPE==1) || (DUMP_TYPE==3) ) {
941 if( (fout = fopen(prefix_dump[hh][count], "w")) == NULL) {
942 fprintf(stderr, "Error opening file %s.",
943 prefix_dump[hh][count]);
944 exit(139);
945 }
946
947 for( bb=0 ; bb<idx2 ; bb++ ){
948 for( rr=0 ; rr<4 ; rr++ )
949 fprintf(fout,"%d\t",temp_arr2[bb][rr]);
950 fprintf(fout,"\n");
951 }
952
953 fclose(fout);
954 }
955 for( bb=0 ; bb<(int) Param_grid[hh][idx3][0] ; bb++) // free all
956 free(temp_arr2[bb]);
957 free(temp_arr2);
958 }
959 }
960 count++;
961 } // end of network loop
962
963 // ****** freeing **********
964 for ( bb = 0 ; bb < N_nets ; bb++ )
965 for ( rr = 0 ; rr < N_totpair[bb] ; rr++ ) {
966 free(prefix_dump[bb][rr]);
967 free(prefix_netmap[bb][rr]);
968 }
969 for ( bb = 0 ; bb < N_nets ; bb++ ) {
970 free(prefix_dump[bb]);
971 free(prefix_netmap[bb]);
972 }
973 free(prefix_dump);
974 free(prefix_netmap);
975
976 free(networkMAPS);
977 }
978
979 free(N_totpair);
980 /*if(str_lab1) {
981 str_lab1 = NULL;
982 free(str_lab1);
983 }
984 if(str_lab2) {
985 str_lab2 = NULL;
986 free(str_lab2);
987 }*/
988
989 RETURN(1);
990 }
991
Setup_Labels_Indices_Unc_M_both(int * Dim,int *** mskd,int *** INDEX,int *** INDEX2,float ** UNC,float ** coorded,float ** copy_coorded,THD_3dim_dataset * insetFA,short * DirPerVox,int N_HAR,THD_3dim_dataset ** insetV,THD_3dim_dataset * insetUC,float unc_minei_std,float unc_minfa_std,int N_nets,int * NROI,THD_3dim_dataset * mset1,int ** MAPROI,int ** INV_LABELS,int *** NETROI)992 int Setup_Labels_Indices_Unc_M_both(int *Dim, int ***mskd, int ***INDEX,
993 int ***INDEX2, float **UNC,
994 float **coorded, float **copy_coorded,
995 THD_3dim_dataset *insetFA,
996 short *DirPerVox,
997 int N_HAR,
998 THD_3dim_dataset **insetV,
999 THD_3dim_dataset *insetUC,
1000 float unc_minei_std, float unc_minfa_std,
1001 int N_nets, int *NROI,
1002 THD_3dim_dataset *mset1, int **MAPROI,
1003 int **INV_LABELS, int ***NETROI)
1004 {
1005
1006 int i,j,k,idx,m,mm,rr;
1007 int idw,n,nn;
1008 float tempvmagn;
1009 float temp1;
1010
1011 // set up eigvecs in 3D coord sys,
1012 // mark off where ROIs are and keep index handy
1013 // also make uncert matrix, with min values of del ang and del FA
1014 for( k=0 ; k<Dim[2] ; k++ )
1015 for( j=0 ; j<Dim[1] ; j++ )
1016 for( i=0 ; i<Dim[0] ; i++ )
1017 if(mskd[i][j][k]) {
1018 idx = INDEX2[i][j][k];
1019 idw = INDEX[i][j][k];
1020
1021 coorded[idx][0]=copy_coorded[idx][0]=
1022 THD_get_voxel(insetFA,idw,0);
1023
1024 // Uncertainty parts if
1025 if( UNC != NULL ) {
1026 if( N_HAR==0 ) { // DTI
1027 // all from data set
1028 //for( m=0 ; m<6 ; m++ )
1029 // UNC[idx][m] = THD_get_voxel(insetUC,idw,m);
1030 temp1 = (THD_get_voxel(insetUC,idw,1) > unc_minei_std ) ?
1031 THD_get_voxel(insetUC,idw,1) : unc_minei_std;
1032 temp1 = pow(temp1,2);
1033 temp1+= pow(THD_get_voxel(insetUC,idw,0),2);
1034 UNC[idx][0] = sqrt(temp1); // delta e_{12}
1035
1036 temp1 = (THD_get_voxel(insetUC,idw,3) > unc_minei_std ) ?
1037 THD_get_voxel(insetUC,idw,3) : unc_minei_std;
1038 temp1 = pow(temp1,2);
1039 temp1+= pow(THD_get_voxel(insetUC,idw,2),2);
1040 UNC[idx][1] = sqrt(temp1); // delta e_{13}
1041
1042 // mean and std of delta FA
1043 UNC[idx][2] = THD_get_voxel(insetUC,idw,4);
1044 UNC[idx][3] = (THD_get_voxel(insetUC,idw,5) > unc_minfa_std ) ?
1045 THD_get_voxel(insetUC,idw,5) : unc_minfa_std;
1046
1047 }
1048 else { // HARDI
1049 // all from data set
1050 for( m=0 ; m<DirPerVox[idx] ; m++ ) {
1051 UNC[idx][m] =
1052 (THD_get_voxel(insetUC,idw,m+1) > unc_minei_std ) ?
1053 THD_get_voxel(insetUC,idw,m+1) : unc_minei_std;
1054 }
1055 UNC[idx][m] =
1056 (THD_get_voxel(insetUC,idw,0) > unc_minfa_std) ?
1057 THD_get_voxel(insetUC,idw,0) : unc_minfa_std;
1058 }
1059
1060 }
1061
1062 // setup all vectors
1063 for( n=0 ; n<DirPerVox[idx] ; n++ ) {
1064 nn = 3*n;
1065
1066 for( m=0 ; m<3 ; m++ )
1067 coorded[idx][nn+m+1]=copy_coorded[idx][nn+m+1]=
1068 THD_get_voxel(insetV[n],idw,m);
1069
1070 // apparently, some |V1| != 1... gotta fix
1071 tempvmagn = sqrt(copy_coorded[idx][nn+1]*copy_coorded[idx][nn+1]+
1072 copy_coorded[idx][nn+2]*copy_coorded[idx][nn+2]+
1073 copy_coorded[idx][nn+3]*copy_coorded[idx][nn+3]);
1074 // if(tempvmagn<0.99)
1075 if(tempvmagn>0)
1076 for( m=1 ; m<4 ; m++ ) {
1077 copy_coorded[idx][nn+m]/=tempvmagn;
1078 coorded[idx][nn+m]/=tempvmagn;
1079 }
1080 }
1081
1082 for( m=0 ; m<N_nets ; m++ ) {
1083 // allow indentification by index --
1084 // now, since we allow non-consecutive ROI labels,
1085 // just use the compacted list numbers via INV_LABELS
1086 if( THD_get_voxel(mset1, idw, m)>0.5 )
1087 MAPROI[idx][m] = INV_LABELS[m][(int) THD_get_voxel(mset1,idw,m)];
1088 else if( THD_get_voxel(mset1, idw, m)<-0.5 )
1089 // silly, is zero anyways... can use this to set to neg mask
1090 MAPROI[idx][m] = -1;
1091
1092 // ppppretty sure there's no need for this here->
1093 // it's all zeros already
1094 // counter for number of kept tracks passing through
1095 //rr = NROI[m]*(NROI[m]+1);
1096 //rr/= 2;
1097 //for( mm=0 ; mm<rr ; mm++ )
1098 // NETROI[idx][m][mm]) = 0;
1099 }
1100 }
1101
1102 RETURN(1);
1103
1104 }
1105
1106
Setup_Ndir_per_vox(int N_HAR,int * Dim,int *** mskd,int *** INDEX,int *** INDEX2,THD_3dim_dataset ** insetHARDIR,short * DirPerVox)1107 int Setup_Ndir_per_vox(int N_HAR, int *Dim, int ***mskd,
1108 int ***INDEX,
1109 int ***INDEX2,
1110 THD_3dim_dataset **insetHARDIR,
1111 short *DirPerVox)
1112 {
1113 int i,j,k,idx,m,n;
1114 float tempvmagn;
1115
1116 // set up eigvecs in 3D coord sys,
1117 // mark off where ROIs are and keep index handy
1118 // also make uncert matrix, with min values of del ang and del FA
1119 for( k=0 ; k<Dim[2] ; k++ )
1120 for( j=0 ; j<Dim[1] ; j++ )
1121 for( i=0 ; i<Dim[0] ; i++ )
1122 if(mskd[i][j][k]) {
1123 idx = INDEX2[i][j][k];
1124 for( n=0 ; n<N_HAR ; n++ ){
1125 tempvmagn = 0;
1126 for( m=0 ; m<3 ; m++ )
1127 tempvmagn+=
1128 THD_get_voxel(insetHARDIR[n],INDEX[i][j][k],m)*
1129 THD_get_voxel(insetHARDIR[n],INDEX[i][j][k],m);
1130
1131 if( tempvmagn > 0.01)
1132 DirPerVox[idx]+=1;
1133 else if( tempvmagn <0.00001 ) // i.e., zero
1134 continue;
1135 else { // awkwardly small.
1136 INFO_message("Minor note: there is a tiny"
1137 " (magn < 0.1)"
1138 " vector in the %d-th direction set."
1139 "\n\t--> Will exclude that voxel"
1140 " from walkable mask-- recommend checking"
1141 " model fit.",n);
1142 mskd[i][j][k] = 0;
1143 }
1144 }
1145 }
1146 return 1;
1147 }
1148
1149 /*
1150 For Monte Carlo iterations, perturb values in mask
1151 */
DTI_Perturb_M(int * Dim,int *** mskd,int *** INDEX,int *** INDEX2,float ** UNC,float ** coorded,float ** copy_coorded,gsl_rng * r,THD_3dim_dataset ** insetV)1152 int DTI_Perturb_M( int *Dim, int ***mskd, int ***INDEX, int ***INDEX2,
1153 float **UNC, float **coorded, float **copy_coorded,
1154 gsl_rng *r,
1155 THD_3dim_dataset **insetV)
1156 {
1157
1158 int i,j,k,m,idx,idw;
1159 float tempvmagn,testang,w2,w3;
1160 float tempv[3];
1161
1162 // have already got `del theta' for angles in UNC!
1163
1164 for( k=0 ; k<Dim[2] ; k++ )
1165 for( j=0 ; j<Dim[1] ; j++ )
1166 for( i=0 ; i<Dim[0] ; i++ ) {
1167 idx = INDEX2[i][j][k];
1168 idw = INDEX[i][j][k];
1169 // only do region in brain mask
1170 if(mskd[i][j][k]) {
1171
1172 // these are weights determined by rotation angle,
1173 // (prob. determined by jackknifing with 3dDWUncert)
1174 // each tips in the +/- direc toward/away from each evec
1175 // by averaging and that's why tan of angle is taken
1176 //@@@
1177 testang = GSL_RAN(r,1.0)*UNC[idx][0];
1178 w2 = tan(testang);
1179
1180 testang = GSL_RAN(r,1.0)*UNC[idx][1];
1181 w3 = tan(testang);
1182
1183
1184 for( m=0 ; m<3 ; m++)
1185 tempv[m] = coorded[idx][m+1] +
1186 w2*THD_get_voxel(insetV[1],idw,m) +
1187 w3*THD_get_voxel(insetV[2],idw,m);
1188 tempvmagn = sqrt(tempv[0]*tempv[0]+
1189 tempv[1]*tempv[1]+tempv[2]*tempv[2]);
1190
1191 for( m=0 ; m<3 ; m++)
1192 copy_coorded[idx][m+1] = tempv[m]/tempvmagn;
1193
1194 // apply bias and std of FA
1195 copy_coorded[idx][0] = coorded[idx][0] + UNC[idx][2] +
1196 ( UNC[idx][3] * GSL_RAN(r,1.0) );
1197
1198 }
1199 }
1200
1201 RETURN(1);
1202 }
1203
1204
HARDI_Perturb(int * Dim,int *** mskd,int *** INDEX,int *** INDEX2,float ** UNC,float ** coorded,float ** copy_coorded,gsl_rng * r,short * DirPerVox)1205 int HARDI_Perturb( int *Dim, int ***mskd, int ***INDEX, int ***INDEX2,
1206 float **UNC, float **coorded, float **copy_coorded,
1207 gsl_rng *r, short *DirPerVox)
1208 {
1209
1210 int i,j,k,m,idx,n,B,nn;
1211 float tempvmagn,testang,w2,w3;
1212 float for_pol, for_azim;
1213 float tempv[3];
1214 float rot[3][3];
1215
1216 for( k=0 ; k<Dim[2] ; k++ )
1217 for( j=0 ; j<Dim[1] ; j++ )
1218 for( i=0 ; i<Dim[0] ; i++ ) {
1219 idx = INDEX2[i][j][k];
1220 // only do region in brain mask
1221 if(mskd[i][j][k]) {
1222
1223 for( n=0 ; n<DirPerVox[idx] ; n++ ){
1224 nn=3*n+1; // plus one b/c zeroth brick is FA
1225
1226 // unit vect, randomly perturbed around (0,0,1)
1227 for_pol = GSL_RAN(r,1.0)*UNC[idx][n];
1228 for_azim = TWOPI*gsl_rng_uniform (r); // rand in 2*pi*[0,1)
1229 tempv[0] = tempv[1] = sin(for_pol);
1230 tempv[0]*= cos(for_azim);
1231 tempv[1]*= sin(for_azim);
1232 tempv[2] = cos(for_pol);
1233
1234 for_pol = acos( coorded[idx][nn+2] ); // acos(z)
1235 for_azim = atan2( coorded[idx][nn+1],coorded[idx][nn+0]); // atan(y,x)
1236
1237 // do rotation: select within func for which vec
1238 // (nn) because 0th brick is now FA value
1239 B = Two_DOF_Rot(tempv, copy_coorded[idx]+nn,
1240 for_pol, for_azim, rot);
1241
1242 /* printf("\n%d DP = %f; for (%f, %f, %f) and (%f, %f, %f)",n,
1243 copy_coorded[idx][nn+0]*coorded[idx][nn+0]+
1244 copy_coorded[idx][nn+1]*coorded[idx][nn+1]+
1245 copy_coorded[idx][nn+2]*coorded[idx][nn+2],
1246 coorded[idx][nn+0],
1247 coorded[idx][nn+1],
1248 coorded[idx][nn+2],
1249 copy_coorded[idx][nn+0],
1250 copy_coorded[idx][nn+1],
1251 copy_coorded[idx][nn+2]
1252 );*/
1253 }
1254
1255 // for FA; n should have correct value here...
1256 //for_pol = GSL_RAN(r,1.0)*UNC[idx][n];
1257 copy_coorded[idx][0] = coorded[idx][0] +
1258 GSL_RAN(r,1.0)*UNC[idx][n];
1259 }
1260 }
1261
1262 RETURN(1);
1263 }
1264
1265
Two_DOF_Rot(float * X,float * Y,double POL,double AZIM,float rot[3][3])1266 int Two_DOF_Rot(float *X, float *Y,
1267 double POL, double AZIM, float rot[3][3] )
1268 {
1269 int i, j, k;
1270 float C0, C1, S0, S1;
1271
1272 C0 = cos(POL);
1273 S0 = sin(POL);
1274 C1 = cos(AZIM);
1275 S1 = sin(AZIM);
1276
1277 rot[0][0] = C0*C1; rot[0][1] = -S1; rot[0][2] = C1*S0;
1278 rot[1][0] = C0*S1; rot[1][1] = C1; rot[1][2] = S0*S1;
1279 rot[2][0] = -S0; rot[2][1] = 0.; rot[2][2] = C0;
1280
1281 for( i=0 ; i<3 ; i++)
1282 Y[i] = 0;
1283
1284 for( i=0 ; i<3 ; i++)
1285 for( j=0 ; j<3 ; j++)
1286 Y[i]+= rot[i][j]*X[j];
1287
1288 RETURN(1);
1289 }
1290
1291
1292 /*
1293 NEW NEW ONE. Multi directionality.
1294 */
1295
TrackItP_NEW_M(int NHAR,short * DirPerVox,int SEL,float ** CC,int * IND,float * PHYSIND,float * Edge,int * dim,float minFA,float maxAng,int arrMax,int ** T,float ** flT,int FB,float * physL,int *** ID2)1296 int TrackItP_NEW_M(int NHAR, short *DirPerVox, int SEL, float **CC,
1297 int *IND, float *PHYSIND,
1298 float *Edge, int *dim, float minFA,
1299 float maxAng, int arrMax, int **T,
1300 float **flT, int FB, float *physL,
1301 int ***ID2)
1302 {
1303 int tracL = 0; // trac length
1304 float Iam0[3]; // 'physical' location
1305 float AA, BB; // mult by L: (1-f)/2 and (1+f)/2, resp., define diagonal
1306 int vsign[3]; // keep trac of vel component sign
1307 float test[3]; // use to see if we move diag or not
1308 int ord[3];
1309 float dotprod;
1310 int go[3]; // for walking through
1311 int win;
1312 float targedge[3]; // possible walls to aim for
1313 float stest[3]; // compare s values-- get shortest 'time' to wall
1314 int m,n,tt;
1315 int walkback;
1316 float physdist = 0.0; // init dist walked is 0;
1317 float FF = 0.4143; //from: (sin(22.5*CONV)/sin((90-22.5)*CONV));
1318 float divid;
1319 int FULLSEL;
1320 float which_dp;
1321
1322 ENTRY("TrackItP_NEW_M");
1323
1324 AA = 0.5*(1.0-FF);
1325 BB = 0.5*(1.0+FF);
1326
1327 // initial place in center of first volume
1328 for( n=0 ; n<3 ; n++)
1329 Iam0[n] = PHYSIND[n];
1330 tt = ID2[IND[0]][IND[1]][IND[2]]; // tt changes when INDs do
1331 // init dotprod is ~unity
1332 dotprod = 0.999;
1333 // square we are 'in' which survived tests: keep trac of, and add
1334 // to temp list
1335 for( n=0 ; n<3 ; n++) {
1336 T[tracL][n] = IND[n];
1337 flT[tracL][n] = Iam0[n];
1338 }
1339 tracL+= 1;
1340
1341 // tracking!
1342 // conditions to stop are:
1343 // + too long of a tract (most likely something wrong with alg,
1344 // because max is large)
1345 // + FA drops below threshold
1346 // + angulation goes above threshold
1347 while( (tracL < arrMax) && (CC[tt][0] >= minFA)
1348 && ( dotprod >= maxAng) ) {
1349
1350 // offset in arrays to select which vect because of having:
1351 // 1) FA map as 0th brick in CC,
1352 // 2) possibly hardi/multi
1353 FULLSEL = 3*SEL+1;
1354
1355 // go to nearest edge
1356 for( n=0 ; n<3 ; n++) {
1357 go[n] = 0; // just resetting our direction to
1358 // Designate up/down, L/R, forw/back, with FB value given before.
1359 // Use SEL->FULLSEL to select which vect
1360 if(CC[tt][FULLSEL+n]*FB >=0) {
1361 targedge[n] = (IND[n]+1)*Edge[n]; // physical units
1362 vsign[n] = 1;
1363 }
1364 else {
1365 targedge[n] = IND[n]*Edge[n];
1366 vsign[n] = -1;
1367 }
1368 }
1369
1370
1371 // calc 'param' to get to edge...
1372 for( n=0 ; n<3 ; n++) {
1373 if( fabs(CC[tt][FULLSEL+n]) < EPS_V)
1374 divid = EPS_V*vsign[n];
1375 else
1376 divid = FB*CC[tt][FULLSEL+n];
1377 stest[n] = (targedge[n]-Iam0[n])/divid;
1378 }
1379 walkback=0;
1380
1381 // say, due to very small CC in opp direction, or trying to push
1382 // us back into previous
1383 for( n=0 ; n<3 ; n++)
1384 if( (stest[n]<0) )
1385 walkback =1;
1386
1387
1388 if(walkback==0) {
1389 for( n=0 ; n<3 ; n++) // try this config as initial guess
1390 ord[n] = n;
1391
1392 if(stest[ ord[1] ]<stest[ ord[0] ]) { // switch
1393 ord[0] = 1; // these are known values of each...
1394 ord[1] = 0;
1395 }
1396 if(stest[ ord[2] ]<stest[ ord[0] ]) { // switch
1397 n = ord[0]; // save temp
1398 ord[0] = ord[2]; // overwrite
1399 ord[2] = n; // finish switch
1400 }
1401 if(stest[ ord[2] ]<stest[ ord[1] ]) { // switch
1402 n = ord[1];
1403 ord[1] = ord[2];
1404 ord[2] = n;
1405 }
1406
1407 win = ord[0];
1408 go[ord[0]] = vsign[ord[0]];
1409
1410 // winner is here; other 2 indices haven't changed, test them.
1411 test[ord[1]] = Iam0[ord[1]] +
1412 stest[ord[0]]*FB*CC[tt][FULLSEL+ord[1]] -
1413 (IND[ord[1]]*Edge[ord[1]]);
1414
1415 if( ( (vsign[ord[1]]>0 ) && (test[ord[1]] > Edge[ord[1]]*BB) ) ||
1416 ( (vsign[ord[1]]<0 ) && (test[ord[1]] < Edge[ord[1]]*AA) ) ){
1417 // then test and see where it would end up
1418 test[ord[0]] = Iam0[ord[0]] +
1419 stest[ord[1]]*FB*CC[tt][FULLSEL+ord[0]] -
1420 (IND[ ord[0]] + go[ord[0]])*Edge[ord[0]];
1421
1422 if( ( (vsign[ord[0]]>0) && (test[ord[0]] < Edge[ord[0]]*AA) ) ||
1423 ( (vsign[ord[0]]<0) && (test[ord[0]] > Edge[ord[0]]*BB) ) ){
1424 go[ord[1]] = vsign[ord[1]]; // partially-'diagonal' route
1425 win = ord[1];
1426
1427 // and only now, do we test for the other diagonal
1428 test[ord[2]] = Iam0[ord[2]] +
1429 stest[ord[0]]*FB*CC[tt][FULLSEL+ord[2]] -
1430 (IND[ord[2]]*Edge[ord[2]]);
1431
1432 if(((vsign[ord[2]]>0 ) && (test[ord[2]] > Edge[ord[2]]*BB)) ||
1433 ((vsign[ord[2]]<0 ) && (test[ord[2]] < Edge[ord[2]]*AA)) ){
1434 test[ord[0]] = Iam0[ord[0]] +
1435 stest[ord[2]]*FB*CC[tt][FULLSEL+ord[0]] -
1436 (IND[ord[0]]+go[ord[0]])*Edge[ord[0]];
1437 test[ord[1]] = Iam0[ord[1]] +
1438 stest[ord[2]]*FB*CC[tt][FULLSEL+ord[1]] -
1439 (IND[ord[1]] + go[ord[1]])*Edge[ord[1]];
1440
1441 // check both for diag-diag
1442 if(((vsign[ord[0]]>0) && (test[ord[0]] < Edge[ord[0]]*AA)) ||
1443 ((vsign[ord[0]]<0) && (test[ord[0]] > Edge[ord[0]]*BB)))
1444 if(((vsign[ord[1]]>0) &&
1445 (test[ord[1]] < Edge[ord[1]]*AA)) ||
1446 ((vsign[ord[1]]<0) &&
1447 (test[ord[1]] > Edge[ord[1]]*BB))){
1448 go[ord[2]] = vsign[ord[2]]; // fully-'diagonal' route
1449 win = ord[2];
1450 }
1451 }
1452 }
1453 }
1454
1455 // move to boundary of next square, updating square we are 'in'
1456 // with current eigenvec
1457 for( n=0 ; n<3 ; n++) // phys loc
1458 Iam0[n]+= stest[ win ]*FB*CC[tt][FULLSEL+n];
1459 for( n=0 ; n<3 ; n++) // update indices of square we're in
1460 IND[n] = IND[n]+go[n]; // will change after testing vals
1461
1462 physdist+= stest[win];
1463
1464 // one way we can stop is by trying to 'walk out' of the volume;
1465 // can check that here
1466 if((IND[0] < dim[0]) && (IND[1] < dim[1]) && (IND[2] < dim[2]) &&
1467 (IND[0] >= 0) && (IND[1] >= 0) && (IND[2] >= 0) ) {
1468 tt = ID2[IND[0]][IND[1]][IND[2]];
1469
1470 // dot prod for stopping cond (abs value)
1471 // check with current vec dotproducted with previous
1472 // ++ NOW also checking which vector to select if there are
1473 // multi directions: largest |dotprod| means smallest angle
1474 which_dp = 0.;
1475 SEL = 0;
1476 for( m=0 ; m<DirPerVox[tt] ; m++) {
1477 dotprod = 0.;
1478 for( n=0 ; n<3 ; n++)
1479 dotprod+= CC[tt][3*m+1+n]*
1480 FB*CC[ID2[IND[0]-go[0]][IND[1]-go[1]][IND[2]-go[2]]][FULLSEL+n];
1481 if (fabs(dotprod) > which_dp)
1482 which_dp = fabs(dotprod) ;
1483 SEL = m;
1484 }
1485
1486 // because of ambiguity of direc/orient of evecs
1487 // and will be checked for criterion at start of next while loop
1488 // because we will keep moving in 'negative' orientation of evec
1489 if(dotprod<0) {
1490 dotprod*=-1.;
1491 FB = -1;
1492 }
1493 else
1494 FB = 1; // move along current orientation of next one
1495
1496 // make sure we haven't been here before
1497 for( n=0 ; n<tracL ; n++)
1498 if( (IND[0]==T[n][0]) && (IND[1]==T[n][1]) && (IND[2]==T[n][2]) )
1499 dotprod =-1.;
1500
1501 if(dotprod<0) // bad
1502 for( n=0 ; n<3 ; n++) {
1503 T[tracL][n] = -1;
1504 flT[tracL][n] = -1;
1505 }
1506 else // fine
1507 for( n=0 ; n<3 ; n++) {
1508 T[tracL][n] = IND[n];
1509 flT[tracL][n] = Iam0[n];
1510 }
1511 tracL+= 1;
1512
1513 }
1514 else {
1515 // to not try to access inaccessible value
1516 // so we will exit tracking in this direction
1517 // at start of next loop
1518 for( n=0 ; n<3 ; n++) {
1519 IND[n] = 0;
1520 T[tracL][n] = -1;
1521 flT[tracL][n] = -1;
1522 }
1523 dotprod = -2.;
1524 tt=0;
1525 tracL+= 1;
1526 }
1527 }
1528 else{
1529 dotprod = -3.;
1530 tt=0;
1531 for( n=0 ; n<3 ; n++) {
1532 T[tracL][n] = -1;
1533 flT[tracL][n] = -1;
1534 }
1535 tracL+= 1;
1536 }
1537
1538 }
1539
1540 if(dotprod>=0 && (CC[tt][0] >= minFA)){
1541 for( n=0 ; n<3 ; n++) {
1542 T[tracL][n] = -1;
1543 flT[tracL][n] = -1;
1544 }
1545 tracL+= 1;
1546 }
1547 if(tracL >= arrMax) //{
1548 WARNING_message("Err in data set or run conditions:\n"
1549 "\t half-tract reached max arr len! Truncating it.\n");
1550 //exit(101);
1551 // }
1552
1553 physL[0] = physdist;
1554
1555 RETURN(tracL);
1556 }
1557
1558
1559