1 /*
2 Probabilistic(+deterministic) tracking, first draft at
3 AFNIfication, March 2012. using FACTID code, from Taylor,
4 Kuan-Hung, Lin and Biswal (2012)
5
6 Sept 2012: fixing some memory stuff.
7
8 Nov 2012: ability to use ROI labels which are non-consecutive ints.
9
10 Dec 2012:
11 + allow non-FA map to define `WM' for tracks,
12 + can include brainmask, if wanted.
13 + allow both nifti and afni inputs to be read
14
15 Jan 2013:
16 + 'NOT' masks
17 + minimum uncert
18
19 Jan 2013:
20 + output ascii format 'dumps' to be turned into masks by 3dUndump,
21 or individ brain mask dumps, or both
22
23 Feb 2013:
24 + allow for output of unthresholded values for ROI maps
25
26 Mar 2013:
27 + output *.grid file now has row of ROI labels in it
28 + floor instead of ceil to calc NmNsThr
29
30 Aug 2013:
31 + merge with det tracking to have networks here, and single prog: @)
32 + mini prob tracking for det case
33 + better visualization tools, a la Ziad Saad
34 + groundwork for adding SUMA+matrix capabilities, new IO: @)))
35 + will start voxelwise track info
36
37 Sep 2013:
38 + include better NOT mask behavior-- simpler; instead of separate mask,
39 have negative (-1) ROIs in same brik with positive ones
40
41 Jan 2014:
42 + firstly, migrated HERE.
43 + multidirectional tracking: HARDI!
44 + nb: *coorded* now have 0th brick as FA info, and then vec1, vec2, etc.
45 + lots of interior changes: searching for names/globbing, allow more inputs
46 + allow commandline opts for algopts; got rid of M and bval in algopt- unused
47 + change output numbering to start from 0.
48 + Instatract stuff better, ZSS work.
49 + '-cut_at_rois' now default behavior; added '-uncut_at_rois' switch
50 + can switch off *.trk format
51
52 Feb 2014:
53 + couple option name switches, and how TROSS name is saved
54
55 Mar 2014:
56 + '-dti_list' input option added at user request.
57
58 May/June 2014:
59 + structified call to tracking, for ZSS
60 + redoing what reaching arrMax does: throw an error, but don't stop
61 + more stats in grid: frac and phys volumes included
62 + extra through mask ('-thru_mask')
63 + some nicer output-- file name prefix at end, count miniprob
64
65 Aug 2014:
66 + add option for output of PAIRMAP to not be 2^X for connections
67 + add option for some scaling of tracks by ROI vols and surf areas
68
69 Aug 2014 (again):
70 + using UHTs to store information-- major mem savings
71 + prob a bit faster now, as well. No change in results
72
73 Sept 2014:
74 + use label table of netrois, if exists
75
76 Jan 2015:
77 + '-nifti' switch for NIFTI output
78
79 Feb 2015:
80 + option to threshold minip results by bundle size
81 + trk output no longer default
82 + *neeed* logic in DET|MINIP modes
83
84 Oct/Nov 2016:
85 + major streamlining of internal code (DTI globbing)
86 + newer functions to clean up functionality
87
88 Nov 2016:
89 + new distance calculation: output lengths of bundles
90 + new targ_surf_stop opt: stop tracts at surface of nets
91 + new targ_surf_twist opt: stop tracts *before*
92 surface of nets
93 + *lots* of new internal streamlining as part
94 of this process
95
96 Nov 2016b:
97 + got rid of occasional nan in sBL (rounding err neg -> sqrt ->
98 nan) with IF
99
100 Aug 2017:
101 + change readglob function to get only 3D vol files where
102 appropriate-- previous condition would get 1D files and give
103 error
104
105 June 2019:
106
107 + new opt -trk_opp_orient: for TRK file output, one can oppositize
108 the voxel_order setting. This is for using TRK files in other
109 software, such as DTK.
110
111 [PT: July 15, 2020] from include "suma_suma.h" -> "suma_objs.h"
112
113 */
114
115
116 #include <stdio.h>
117 #include <stdlib.h>
118 #include <string.h>
119 #include <math.h>
120 #include <unistd.h>
121 #include <time.h>
122 #include <debugtrace.h>
123 #include "mrilib.h"
124 #include "3ddata.h"
125 #include <gsl/gsl_randist.h>
126 #include <gsl/gsl_rng.h>
127 #include "Fat_Labels.h"
128 #include "DoTrackit.h"
129 #include "TrackIO.h"
130 #include "FuncTrac.h"
131 #include "readglob.h"
132 #include "checks_and_balances.h"
133
134 #include "suma_objs.h"
135
136
137 //#define MAX_PARAMS (11) // for *.trk output
138 //#define N_dti_scal (3)
139 //#define N_dti_scal_tot (6) // 1+2+3, matches with N_dti_scal
140 //#define N_dti_vect (3)
141 //#define N_plus_dtifile (4) // num of allowed additional param files in NIMLver
142 //#define PreCalcParLabs (5) // we define a few things initially
143
144 #define N_DEF_PAR_LABS (7) // we define a few things initially
145 #define N_XTR_PAR_LABS (3) // we define a few things initially
146
147 static char *DEF_PAR_LABS[N_DEF_PAR_LABS] = { "NT", "fNT", "PV",
148 "fNV", "NV", "BL",
149 "sBL"};
150 static char *XTR_PAR_LABS[N_XTR_PAR_LABS] = { "NTpTarVol",
151 "NTpTarSA",
152 "NTpTarSAFA"};
153
154
155
156 int GetSizesOfTargets( int ***ROI_SIZES,
157 int **MAPROI, int *Dim,
158 int ***mskd,
159 int ***INDEX2,
160 float **coorded, float minFA,
161 int N_nets, int *NROI,
162 int **INV_LABELS);
163
164
165 int RunTrackingMaestro( int comline, TRACK_RUN_PARAMS opts,
166 int argc, char *argv[]);
167
168
169
usage_TrackID(int detail)170 void usage_TrackID(int detail)
171 {
172 printf("\n"
173 " FACTID-based tractography code, from Taylor, Cho, Lin and Biswal (2012),\n"
174 " and part of FATCAT (Taylor & Saad, 2013) in AFNI. Version 2.1 (Jan. 2014),\n"
175 " written by PA Taylor and ZS Saad.\n"
176 "\n"
177 " Estimate locations of WM associated with target ROIs, particularly between\n"
178 " pairs of GM in a network; can process several networks in a given run.\n"
179 "\n"
180 " Now does both single tract propagation per voxel (as per DTI) and \n"
181 " multi-directional tracking (as in HARDI-type models). Many extra files can\n"
182 " be loaded in for getting quantitative stats in WM-ROIs, mostly done via\n"
183 " search from entered prefixes. Many more switches and options are available\n"
184 " to the user to control the tracking (yay!).\n"
185 " Track display capabilities in SUMA have been boosted and continue to rise\n"
186 " quickly (all courtesy of ZS Saad).\n"
187 "\n"
188 "****************************************************************************\n"
189 "+ NOTE that this program runs in three separate modes, each with its own\n"
190 " subset of commandline options and outputs:\n"
191 " $ 3dTrackID -mode {DET | MINIP | PROB} ... \n"
192 " where DET -> deterministic tracking,\n"
193 " MINIP -> mini-probabilistic tracking,\n"
194 " PROB -> (full) probabilistic tracking.\n"
195 " So, for example, DET and MINIP produce pretty track-image output,\n"
196 " while PROB only provides volumes; MINIP and PROB make use of\n"
197 " tensor uncertainty to produce more robust results than DET; all\n"
198 " produce quantitative statistical output of WM-ROIs; etc. In some cases,\n"
199 " using a combination of all three might even be variously useful in a\n"
200 " particular study.\n"
201 "****************************************************************************\n"
202 " For DTI, this program reads in tensor-related data from, e.g., 3dDWItoDTI,\n"
203 " and also uses results from 3dDWUncert for uncertainty measures when\n"
204 " necessary.\n"
205 "\n"
206 " For HARDI, this program reads in the direction vectors and WM-proxy map \n"
207 " (such as the diffusion anisotropy coefficient, GFA) created by any source-\n"
208 " right now, there's no HARDI modeler in AFNI. Currently known sources which\n"
209 " are reasonably straightforward to use include DSI-Studio (Yeh et al.,\n"
210 " 2010) and Diffusion Toolkit (Wang et al., 2007). An example script of\n"
211 " outputting Qball model data as NIFTI output from the former software is\n"
212 " included in the FATCAT demo set.\n"
213 "\n"
214 " ...And on that note, it is highly recommended for users to check out the\n"
215 " FATCAT demo set, which can be downloaded and unwrapped simply from the\n"
216 " commandline:\n"
217 " $ @Install_FATCAT_Demo\n"
218 " In that demo are data, a number of scripts, and more detailed descriptions\n"
219 " for using 3dTrackID, as well as other programs in the FATCAT litter.\n"
220 " Recommended to always check that one has the most uptodate version.\n"
221 "\n"
222 "****************************************************************************\n"
223 "\n"
224 "+ INPUT NOTES:\n"
225 " NETWORK MAPS, for any '-mode' of track, given as a single- or multi-brik\n"
226 " file via '-netrois':\n"
227 " Each target ROI is defined by the set of voxels with a given integer >0.\n"
228 " Target ROI labels do not have to be purely consecutive.\n"
229 "\n"
230 " Note on vocabulary, dual usage of 'ROI': an (input) network is made up of\n"
231 " *target ROIs*, between/among which one wants to find WM connections; so,\n"
232 " 3dTrackID outputs locations and stats on those calculated *WM-ROIs*.\n"
233 "\n"
234 "****************************************************************************\n"
235 "\n"
236 "+ OUTPUTS, all named using '-prefix INPREF'; somewhat dependent on tracking\n"
237 " mode being utilized ('-mode {DET | MINIP | PROB}').\n"
238 " Because multiple networks can be input simultaneously as a multi-\n"
239 " brik '-netrois ROIS' file, the output prefix will also have a\n"
240 " numerical designation of its network, matching to the brik of\n"
241 " the ROIS file: thus, INPREF_000* goes with ROIS[0], INPREF_001*\n"
242 " with ROIS[1] (if present), etc. This applies with all types of\n"
243 " output files, now described:\n"
244 " 1) *INDIMAP* BRIK files (output in ALL modes).\n"
245 " For each network with N_ROI target ROIs, this is a N_ROI+1 brik file.\n"
246 " 0th brick contains the number of tracts per voxel which passed through\n"
247 " at least one target ROI in that network (and in '-mode PROB', this\n"
248 " number has been thresholded-- see 'alg_Thresh_Frac' below).\n"
249 " If the target ROIs are consecutively labelled from 1 to N_ROI, then:\n"
250 " Each i-th brick (i running from 1 to N_ROI) contains the voxels\n"
251 " through which tracks hitting that i-th target passed; the value of\n"
252 " each voxel is the number of tracks passing through that location.\n"
253 " Else, then:\n"
254 " Each i-th brick contains the voxels through which the tracks\n"
255 " hitting the j-th target passed (where j may or may not equal i; the\n"
256 " value of j is recorded in the brick label: OR_roi_'j'). The target\n"
257 " ROI connectivity is recorded increasing order of 'j'.\n"
258 " For single-ROI inputs (such as a single wholebrain ROI), only the\n"
259 " [0] brick is output (because [1] would be redundant).\n"
260 " 2) *PAIRMAP* BRIK files (output in ALL modes).\n"
261 " (-> This has altered slightly at the end of June, 2014! No longer using\n"
262 " 2^i notation-- made simpler for reading, assuming individual connection\n"
263 " information for calculations was likely obtained more easily with \n"
264 " '-dump_rois {AFNI | BOTH | AFNI_MAP}...)\n"
265 " For each network with N_ROI target ROIs, this is a N_ROI+1 brik file.\n"
266 " 0th brick contains a binary mask of voxels through which passed a\n"
267 " supra-threshold number of tracks (more than 0 for '-mode {DET | MINIP}'\n"
268 " and more than the user-defined threshold for '-mode PROB') between any\n"
269 " pair of target ROIs in that network (by default, these tracks have been\n"
270 " trimmed to only run between ROIs, cutting off parts than dangle outside\n"
271 " of the connection).\n"
272 " If the target ROIs are consecutively labelled from 1 to N_ROI, then:\n"
273 " Each i-th brick (i running from 1 to N_ROI) contains the voxels\n"
274 " through which tracks hitting that i-th target AND any other target\n"
275 " passed; voxels connecting i- and j-th target ROIs have value j, and\n"
276 " the values are summed if a given voxel is in multiple WM ROIs (i.e.,\n"
277 " for a voxel connecting both target ROIs 2 and 1 as well as 2 and 4,\n"
278 " then the value there in brick [2] would be 1 + 4 = 5).\n"
279 " Else, then:\n"
280 " Each i-th brick contains the voxels through which the tracks\n"
281 " hitting the j-th target AND any other target passed (where j may or\n"
282 " may not equal i; the value of j is recorded in the brick label: \n"
283 " AND_roi_'j'). The same voxel labelling and summing rules described\n"
284 " above also apply here.\n"
285 " For single-ROI inputs (such as a single wholebrain ROI), no PAIRMAP\n"
286 " file is output (because it would necessarily be empty).\n"
287 " 3) *.grid ASCII-text file (output in ALL modes).\n"
288 " Simple text file of output stats of WM-ROIs. It outputs the means and\n"
289 " standard deviations of parameter quantities (such as FA, MD, L1, etc.)\n"
290 " as well as counts of tracks and volumes of WM-ROIs. Each matrix is\n"
291 " square, with dimension N_ROI by N_ROI. Like the locations in a standard\n"
292 " correlation matrix, each element reflects associativity with target\n"
293 " ROIs. A value at element (1,3) is the same as that at (3,1) and tells\n"
294 " about the property of a WM-ROI connecting target ROIs 1 and 3 (consider\n"
295 " upper left corner as (1,1)); diagonal elements provide info of tracks\n"
296 " through (at minimum) that single target ROI-- like OR logic connection.\n"
297 " Format of *.grid file is:\n"
298 " Line 1: number of ROIs in network (padded with #-signs)\n"
299 " Line 2: number of output matrices of stats info (padded with #-signs)\n"
300 " Line 3: list of N_ROI labels for that network\n"
301 " Lines following: first line, label of a property (padded with #), and \n"
302 " then N_ROI lines of the N_ROI-by-N_ROI matrix of that\n"
303 " property;\n"
304 " /repeat/\n"
305 " The first *five* matrices are currently (this may change over time):\n"
306 " NT = number of tracks in that WM-ROI\n"
307 " fNT = fractional number of tracks in that WM-ROI, defined as NT\n"
308 " divided by total number of tracts found (may not be relevant)\n"
309 " PV = physical volume of tracks, in mm^3\n"
310 " fNV = fractional volume of tracks compared to masked (internally or\n"
311 " '-mask'edly) total volume; would perhaps be useful if said\n"
312 " mask represents the whole brain volume well.\n"
313 " NV = number of voxels in that WM-ROI.\n"
314 " BL = average length (in mm) of a bundle of tracts.\n"
315 " sBL = stdev of the length (in mm) of a bundle of tracts.\n"
316 " Then, there can be a great variety in the remaining matrices, depending\n"
317 " on whether one is in DTI or HARDI mode and how many scalar parameter\n"
318 " files get input (max is 10). For each scalar file there are two\n"
319 " matrices: first a label (e.g., 'FA') and then an N_ROI-by-N_ROI matrix\n"
320 " of the means of that parameter in each WM-ROI; then a label (here,\n"
321 " would be 'sFA') and then an N_ROI-by-N_ROI matrix of the standard\n"
322 " deviations of that parameter in each WM-ROI.\n"
323 " 4) *niml.tract NIML/SUMA-esque file (output in '-mode {DET | MINIP}')\n"
324 " File for viewing track-like output in SUMA, with, e.g.:\n"
325 " $ suma -tract FILE.niml.tract\n"
326 " 5) *niml.dset NIML/SUMA-esque file (output in '-mode {DET | MINIP}')\n"
327 " File accompanying the *.niml.tract file-- also for use in SUMA, for\n"
328 " including GRID-file like information with the tract info.\n"
329 " $ suma -tract FILE.niml.tract -gdset FILE.niml.dset\n"
330 " 6) *.trk TrackVis-esque file (output in '-mode {DET | MINIP}')\n"
331 " File for viewing track-like output in TrackVis (separate install from\n"
332 " AFNI/SUMA); things mainly done via GUI interface.\n"
333 "\n"
334 "****************************************************************************\n"
335 "\n"
336 "+ LABELTABLE LABELLING (Sept 2014).\n"
337 " The ability to use label tables in tracking result output has been\n"
338 " included. \n"
339 " Default behavior will be to *construct* a labeltable from zero-padded ints\n"
340 " in the '-netrois' file which define target ROIs. Thus, the ROI of '3's\n"
341 " will be given a label '003'. This will be used in INDIMAP and PAIRMAP\n"
342 " brick labels (which is useful if the targets are not consecutively\n"
343 " numbered from 1), PAIRMAP connections in bricks >0, and output \n"
344 " *.niml.tract files. The PAIRMAP labeltable will be created and output\n"
345 " as 'PREFIX_PAIRMAP.niml.lt', and will be useful for the user in (some-\n"
346 " what efficiently) resolving multiple tracts passing through voxels.\n"
347 " These labels are also used in the naming of '-dump_rois AFNI' output.\n"
348 " At the moment, in a given PAIRMAP brick of index >0, labels can only \n"
349 " describe up to two connections through a given voxel. In brick 1, if \n"
350 " voxel is intersected by tracts connection ROIs 1 and 3 as well as ROIs\n"
351 " 1 and 6, then the label there would be '003<->006'; if another voxel\n"
352 " in that brick had those connections as well as one between ROIs 1 and \n"
353 " 4, then the label might be '_M_<->003<->006', or '_M_<->003<->004', or\n"
354 " any two of the connections plus the leading '_M_' that stands for \n"
355 " 'multiple others' (NB: which two are shown is not controlled, but I \n"
356 " figured it was better to show at least some, rather than just the \n"
357 " less informative '_M_' alone). In all of these things, the PAIRMAP\n"
358 " map is a useful, fairly efficient guide-check, but the overlaps are\n"
359 " difficult to represent fully and efficiently, given the possibly\n"
360 " complexity of patterns. For more definite, unique, and scriptable\n"
361 " information of where estimated WM connections are, use the \n"
362 " '-dump_rois AFNI' or '-dump_rois AFNI_MAP' option.\n"
363 " If the '-netrois' input has a labeltable, then this program will program\n"
364 " will read it in, use it in PAIRMAP and INDIMAP bricklabels, PAIRMAP\n"
365 " subbricks with index >0, *niml.tract outputs and, by default, in the\n"
366 " naming of '-dump_rois AFNI' output. The examples and descriptions\n"
367 " directly above still hold, but in cases where the ROI number has an\n"
368 " explicit label, then the former is replaced by the latter's string.\n"
369 " In cases where an input label table does not cover all ROI values, \n"
370 " there is no need to panic-- the explicit input labels will be used\n"
371 " wherever possible, and the zero-padded numbers will be used for the \n"
372 " remaining cases. Thus, one might see PAIRMAP labels such as:\n"
373 " '003<->Right-Amygdala', '_M_<->ctx-lh-insula<->006', etc.\n"
374 "\n"
375 "****************************************************************************\n"
376 "\n"
377 "+ RUNNING AND COMMANDLINE OPTIONS: pick a MODEL and a MODE.\n"
378 " There are now two types of models, DTI and HARDI, that can be tracked.\n"
379 " In HARDI, one may have multiple directions per voxel along which tracts\n"
380 " may propagate; in DTI, there can be only one. Each MODEL has some\n"
381 " required, and some optional, inputs.\n"
382 " Additionally, tracking is run in one of three modes, as described near the\n"
383 " top of this document, '-mode {DET | MINIP | PROB}', for deterministic\n"
384 " mini-probabilistic, or full probabilistic tracking, respectively.\n"
385 " Each MODE has some required, and some optional, inputs. Some options\n"
386 " find work in multiple modes.\n"
387 " To run '3dTrackID', one needs to have both a model and a mode in mind (and\n"
388 " in data...). Below is a table to show the various options available\n"
389 " for the user to perform tracking. The required options for a given\n"
390 " model or mode are marked with a single asterisk (*); the options under\n"
391 " the /ALL/ column are necessary in any mode. Thus, to run deterministic\n"
392 " tracking with DTI data, one *NEEDS* to select, at a minimum:\n"
393 " '-mode DET', '-netrois', '-prefix', '-logic';\n"
394 " and then there is a choice of loading DTI data, with either:\n"
395 " '-dti_in' or '-dti_list',\n"
396 " and then one can also use '-dti_extra', '-mask', '-alg_Nseed_Y',\n"
397 " et al. from the /ALL/ and DET colums; one canNOT specify '-unc_min_FA'\n"
398 " here -> the option is in an unmatched mode column.\n"
399 " Exact usages of each option, plus formats for any arguments, are listed\n"
400 " below. Default values for optional arguments are also described.\n"
401 "\n"
402 " +-----------------------------------------------------------------+\n"
403 " | COMMAND OPTIONS FOR TRACKING MODES AND MODELS |\n"
404 " +-----------------------------------------------------------------+\n"
405 " | /ALL/ | DET | MINIP | PROB |\n"
406 "+--------+-------------------+-------------+-------------+-----------------+\n"
407 " |{dti_in, dti_list}*| | | |\n"
408 " DTI | dti_extra | | | |\n"
409 " | dti_search_NO | | | |\n"
410 "+-~or~---+-------------------+-------------+-------------+-----------------+\n"
411 " | hardi_gfa* | | | |\n"
412 " HARDI | hardi_dirs* | | | |\n"
413 " | hardi_pars | | | |\n"
414 "==~and~==+===================+=============+=============+=================+\n"
415 " | mode* | | | |\n"
416 " OPTIONS | netrois* | | | |\n"
417 " | prefix* | | | |\n"
418 " | mask | | | |\n"
419 " | thru_mask | | | |\n"
420 " | targ_surf_stop | | | |\n"
421 " | targ_surf_twixt | | | |\n"
422 " | | logic* | logic* | |\n"
423 " | | | mini_num* | |\n"
424 " | | | uncert* | uncert* |\n"
425 " | | | unc_min_FA | unc_min_FA |\n"
426 " | | | unc_min_V | unc_min_V |\n"
427 " | algopt | | | |\n"
428 " | alg_Thresh_FA | | | |\n"
429 " | alg_Thresh_ANG | | | |\n"
430 " | alg_Thresh_Len | | | |\n"
431 " | | alg_Nseed_X | alg_Nseed_X | |\n"
432 " | | alg_Nseed_Y | alg_Nseed_Y | |\n"
433 " | | alg_Nseed_Z | alg_Nseed_Z | |\n"
434 " | | | | alg_Thresh_Frac |\n"
435 " | | | | alg_Nseed_Vox |\n"
436 " | | | | alg_Nmonte |\n"
437 " | uncut_at_rois | | | |\n"
438 " | do_trk_out | | | |\n"
439 " | trk_opp_orient | | | |\n"
440 " | dump_rois | | | |\n"
441 " | dump_no_labtab | | | |\n"
442 " | dump_lab_consec | | | |\n"
443 " | posteriori | | | |\n"
444 " | rec_orig | | | |\n"
445 " | tract_out_mode | | | |\n"
446 " | write_opts | | | |\n"
447 " | write_rois | | | |\n"
448 " | pair_out_power | | | |\n"
449 "+--------+-------------------+-------------+-------------+-----------------+\n"
450 "*above, asterisked options are REQUIRED for running the given '-mode'.\n"
451 " With DTI data, one must use either '-dti_in' *or* '-dti_list' for input.\n"
452 "\n"
453 " FOR MODEL DTI:\n"
454 " -dti_in INPREF :basename of DTI volumes output by, e.g., 3dDWItoDT.\n"
455 " NB- following volumes are *required* to be present:\n"
456 " INPREF_FA, INPREF_MD, INPREF_L1,\n"
457 " INPREF_V1, INPREF_V2, INPREF_V3,\n"
458 " and (now) INPREF_RD (**now output by 3dTrackID**).\n"
459 " Additionally, the program will search for all other\n"
460 " scalar (=single brik) files with name INPREF* and will\n"
461 " load these in as additional quantities for WM-ROI\n"
462 " stats; this could be useful if, for example, you have\n"
463 " PD or anatomical measures and want mean/stdev values\n"
464 " in the WM-ROIs (to turn this feature off, see below,\n"
465 " 'dti_search_NO'); all the INPREF* files must be in same\n"
466 " DWI space.\n"
467 " Sidenote: including/omitting a '_' at the end of INPREF\n"
468 " makes no difference in the hunt for files.\n"
469 " -dti_extra SET :if you want to use a non-FA derived definition for the\n"
470 " WM skeleton in which tracts run, you can input one, and\n"
471 " then the threshold in the -algopt file (or, via the\n"
472 " '-alg_Thresh_FA' option) will be applied to \n"
473 " thresholding this SET; similarly for the minimum\n"
474 " uncertainty by default will be set to 0.015 times the\n"
475 " max value of SET, or can be set with '-unc_min_FA'.\n"
476 " If the SET name is formatted as INPREF*, then it will\n"
477 " probably be included twice in stats, but that's not the\n"
478 " worst thing. In grid files, name of this quantity will\n"
479 " be 'XF' (stands for 'extra file').\n"
480 " -dti_search_NO :turn off the feature to search for more scalar (=single\n"
481 " brik) files with INPREF*, for including stats in output\n"
482 " GRID file. Will only go for FA, MD, L1 and RD scalars\n"
483 " with INPREF.\n"
484 " -dti_list FILE :an alternative way to specify DTI input files, where\n"
485 " FILE is a NIML-formatted text file that lists the\n"
486 " explicit/specific files for DTI input. This option is\n"
487 " used in place of '-dti_in' and '-dti_extra' for loading\n"
488 " data sets of FA, MD, L1, etc. An 'extra' set (XF) can\n"
489 " be loaded in the file, as well as supplementary scalar\n"
490 " data sets for extra WM-ROI statistics.\n"
491 " See below for a 'DTI LIST FILE EXAMPLE'.\n"
492 " FOR MODEL HARDI:\n"
493 " -hardi_gfa GFA :single brik data set with generalized FA (GFA) info.\n"
494 " In reality, it doesn't *have* to be a literal GFA, esp.\n"
495 " if you are using some HARDI variety that doesn't have\n"
496 " a specific GFA value-- in such a case, use whatever\n"
497 " could be thresholded as your proxy for WM.\n"
498 " The default threshold is still 0.2, so you will likely\n"
499 " need to set a new one in the '-algopt ALG_FILE' file or\n"
500 " from the commandline with '-alg_Thresh_FA', which does\n"
501 " apply to the GFA in the HARDI case as well.\n"
502 " Stats in GRID file are output under name 'GFA'.\n"
503 " -hardi_dirs DIRS :For tracking if X>1 propagation directions per voxel\n"
504 " are given, for example if HARDI data is input. DIRS\n"
505 " would then be a file with 3*X briks of (x,y,z) ordered,\n"
506 " unit magnitude vector components; i.e., brik [0]\n"
507 " contains V1_x, [1] V1_y, [2] V1_z, [3] V2_x, etc.\n"
508 " (NB: even if X=1, this option works, but that would\n"
509 " seem to take the HAR out of HARDI...)\n"
510 " -hardi_pars PREF :search for scalar (=single brik) files of naming\n"
511 " format PREF*. These will be read in for WM-ROI stats\n"
512 " output in the GRID file. For example, if there are\n"
513 " some files PREF_PD.nii.gz, PREF_CAT.nii.gz and\n"
514 " PREF_DOG.nii.gz, they will be labelled in the GRID file\n"
515 " as 'PD', 'CAT' and 'DOG' (that '_' will be cut out).\n"
516 " MODEL-INDEPENDENT OPTIONS:\n"
517 " -mode MODUS :this necessary option is used to define whether one is\n"
518 " performing deterministic, mini-probabilistic or full-\n"
519 " probabilistic tractography, by selecting one of three\n"
520 " respective modes: DET, MINIP, or PROB.\n"
521 " -netrois ROIS :mask(s) of target ROIs- single file can have multiple\n"
522 " briks, one per network. The target ROIs through which\n"
523 " tracks will be kept should have index values >0. It is\n"
524 " also possible to define anti-targets (exclusionary\n"
525 " regions) which stop a propagating track... in its \n"
526 " tracks. These are defined per network (i.e., per brik)\n"
527 " by voxels with values <0.\n"
528 " -prefix PREFIX :output file name part.\n"
529 " -mask MASK :can include a brainmask within which to calculate \n"
530 " things. Otherwise, data should be masked already.\n"
531 " -thru_mask TM :optional extra restrictor mask, through which paths are\n"
532 " (strictly) required to pass in order to be included\n"
533 " when passing through or connecting targets. It doesn't\n"
534 " discriminate based on target ROI number, so it's\n"
535 " probably mostly useful in examining specific pairwise\n"
536 " connections. It is also not like one of the target\n"
537 " '-netrois' in that no statistics are calculated for it.\n"
538 " Must be same number of briks as '-netrois' set.\n"
539 " -targ_surf_stop :make the final tracts and tracked regions stop at the\n"
540 " outer surface of the target ROIs, rather than being\n"
541 " able to journey arbitrarily far into them (latter being\n"
542 " the default behavior. Might be useful when you want\n"
543 " meaningful distances *between* targets. Tracts stop\n"
544 " after going *into* the outer layer of a target.\n"
545 " This can be applied to DET, MINIP, or PROB modes.\n"
546 " NB: this only affects the connections between pairs\n"
547 " of targets (= AND-logic, off-diagonal elements in\n"
548 " output matrices), not the single-target tracts\n"
549 " (= OR-logic, on-diagonal elements in output\n"
550 " matrices); see also a related option, below.\n"
551 " -targ_surf_twixt :quite similar to '-targ_surf_stop', above, but the\n"
552 " tracts stop *before* entering the target surfaces, so\n"
553 " that they are only between (or betwixt) the targets.\n"
554 " Again, only affects tracts between pairs of targets.\n"
555 "\n"
556 " -logic {OR|AND} :when in one of '-mode {DET | MINIP}', one will look for\n"
557 " either OR- or AND-logic connections among target ROIs\n"
558 " per network (multiple networks can be entered as\n"
559 " separate briks in '-netrois ROIS'): i.e., one keeps\n"
560 " either any track going through at least one network ROI\n"
561 " or only those tracks which join a pair of ROIs.\n"
562 " When using AND here, default behavior is to only keep\n"
563 " voxels in tracks between the ROIs they connect (i.e.,\n"
564 " cut off track bits which run beyond ROIs).\n"
565 " -mini_num NUM :will run a small number NUM of whole brain Monte Carlo\n"
566 " iterations perturbing relevant tensor values in accord\n"
567 " with their uncertainty values (hence, the need for also\n"
568 " using `-uncert' with this option). This might be useful\n"
569 " for giving a flavor of a broader range of connections\n"
570 " while still seeing estimated tracks themselves. NB: if\n"
571 " NUM is large, you might be *big* output track files;\n"
572 " e.g., perhaps try NUM = 5 or 9 or so to start.\n"
573 " Requires '-mode MINIP' in commandline.\n"
574 " -bundle_thr V :the number of tracts for a given connection is called\n"
575 " a bundle. For '-mode {DET | MINIP}', one can choose to\n"
576 " NOT output tracts, matrix info, etc. for any bundle\n"
577 " with fewer than V tracts. This might be useful to weed\n"
578 " out ugly/false tracts (default: V=1).\n"
579 " -uncert U_FILE :when in one of '-mode {MINIP | PROB}', uncertainty\n"
580 " values for eigenvector and WM skeleton (FA, GFA, etc.)\n"
581 " maps are necessary.\n"
582 " When using DTI ('-dti_*'), then use the 6-brik file\n"
583 " from 3dDWUncert; format of the file given below.\n"
584 " When using HARDI ('-hardi_*') with up to X directions\n"
585 " per voxel, one needs U_FILE to have X+1 briks, where\n"
586 " U_FILE[0] is the uncertainty for the GFAfile, and the\n"
587 " other briks are ordered for directions given with\n"
588 " '-hardi_dirs'.\n"
589 " Whatever the values in the U_FILE, this program asserts\n"
590 " a minimum uncertainty of stdevs, with defaults:\n"
591 " for FA it is 0.015, and for GFA or -dti_extra sets it\n"
592 " is 0.015 times the max value present (set with option\n"
593 " '-unc_min_FA');\n"
594 " for each eigenvector or dir, it is 0.06rad (~3.4deg)\n"
595 " (set with option '-unc_min_V')\n"
596 " -unc_min_FA VAL1 :when using '-uncert', one can control the minimum\n"
597 " stdev for perturbing FA (in '-dti_in'), or the EXTRA-\n"
598 " file also in DTI ('-dti_extra'), or GFA (in '-hardi_*).\n"
599 " Default value is: 0.015 for FA, and 0.015 times the max\n"
600 " value in the EXTRA-file or in the GFA file.\n"
601 " -unc_min_V VAL2 :when using '-uncert', one can control the minimum\n"
602 " stdev for perturbing eigen-/direction-vectors.\n"
603 " In DTI, this is for tipping e_1 separately toward e2\n"
604 " and e3, and in HARDI, this is for defining a single\n"
605 " degree of freedom uncertainty cone. Default values are\n"
606 " 0.06rad (~3.4deg) for any eigenvector/direction. User\n"
607 " assigns values in degrees.\n"
608 "\n"
609 " -algopt A_FILE :simple ASCII file with six numbers defining tracking \n"
610 " parameter quantities (see list below); note the\n"
611 " differences whether running in '-mode {DET | MINIP}'\n"
612 " or in '-mode PROB': the first three parameters in each\n"
613 " mode are the same, but the next three differ.\n"
614 " The file can be in the more understandable html-type\n"
615 " format with labels per quantity, or just as a column\n"
616 " of the numbers, necessarily in the correct order.\n"
617 " NB: each quantity can also be changed individually\n"
618 " using a commandline option (see immediately following).\n"
619 " If A_FILE ends with '.niml.opts' (such as would be\n"
620 " produced using the '-write_opts' option), then it is\n"
621 " expected that it is in nice labelled NIML format;\n"
622 " otherwise, the file should just be a column of numbers\n"
623 " in the right order. Examples of A_FILEs are given at\n"
624 " the end of the option section.\n"
625 " -alg_Thresh_FA A :set threshold for DTI FA map, '-dti_extra' FILE, or \n"
626 " HARDI GFA map (default = 0.2).\n"
627 " -alg_Thresh_ANG B :set max angle (in deg) for turning when going to a new\n"
628 " voxel during propagation (default = 60).\n"
629 " -alg_Thresh_Len C :min physical length (in mm) of tracts to keep\n"
630 " (default = 20).\n"
631 " -alg_Nseed_X D :Number of seeds per vox in x-direc (default = 2).\n"
632 " -alg_Nseed_Y E :Number of seeds per vox in y-direc (default = 2).\n"
633 " -alg_Nseed_Z F :Number of seeds per vox in z-direc (default = 2).\n"
634 " +-------> NB: in summation, for example the alg_Nseed_* options\n"
635 " for '-mode {DET | MINIP} place 2x2x2=8 seed points,\n"
636 " equally spread in a 3D cube, in each voxel when\n"
637 " tracking.\n"
638 " -alg_Thresh_Frac G :value for thresholding how many tracks must pass\n"
639 " through a voxel for a given connection before it is\n"
640 " included in the final WM-ROI of that connection.\n"
641 " It is a decimal value <=1, which will multiply the\n"
642 " number of 'starting seeds' per voxel, Nseed_Vox*Nmonte\n"
643 " (see just below for those). (efault = 0.001; for higher\n"
644 " specificity, a value of 0.01-0.05 could be used).\n"
645 " -alg_Nseed_Vox H :number of seeds per voxel per Monte Carlo iteration;\n"
646 " seeds will be placed randomly (default = 5).\n"
647 " -alg_Nmonte I :number of Monte Carlo iterations (default = 1000).\n"
648 " +-------> NB: in summation, the preceding three options for the\n"
649 " '-mode PROB' will mean that 'I' Monte Carlo\n"
650 " iterations will be run, each time using 'H' track\n"
651 " seeds per relevant voxel, and that a voxel will\n"
652 " need 'G*H*I' tracks of a given connection through\n"
653 " it to be included in that WM-ROI. Default example:\n"
654 " 1000 iterations with 5 seeds/voxel, and therefore\n"
655 " a candidate voxel needs at least 0.001*5*1000 = 5\n"
656 " tracks/connection.\n"
657 "\n"
658 " -extra_tr_par :run three extra track parameter scalings for each\n"
659 " target pair, output in the *.grid file. The NT value\n"
660 " of each connection is scaled in the following manners\n"
661 " for each subsequent matrix label:\n"
662 " NTpTarVol :div. by average target volume;\n"
663 " NTpTarSA :div. by average target surface area;\n"
664 " NTpTarSAFA :div. by average target surface area\n"
665 " bordering suprathreshold FA (or equi-\n"
666 " valent WM proxy definition).\n"
667 " NB: the volume and surface area numbers are given in\n"
668 " terms of voxel counts and not using physical units\n"
669 " (consistent: NT values themselves are just numbers.)\n"
670 " -uncut_at_rois :when looking for pairwise connections, keep entire\n"
671 " length of any track passing through multiple targets,\n"
672 " even when part ~overshoots a target (i.e., it's not\n"
673 " between them). When using OR tracking, this is\n"
674 " automatically applied. For probabilistic tracking, not\n"
675 " recommended to use (are untrimmed ends meaningful\?).\n"
676 " The default behavior is to trim the tracts that AND-\n"
677 " wise connect targets to only include sections that are\n"
678 " between the targets, and not parts that run beyond one.\n"
679 " (Not sure why one would want to use this option, to be\n"
680 " honest; see '-targ_surf_stop' for really useful tract\n"
681 " control.)\n"
682 " -dump_rois TYPE :can output individual masks of ROI connections.\n"
683 " Options for TYPE are: {DUMP | AFNI | BOTH | AFNI_MAP}.\n"
684 " Using DUMP gives a set of 4-column ASCII files, each \n"
685 " formatted like a 3dmaskdump data set; it can be recon-\n"
686 " stituted using 3dUndump. Using AFNI gives a set of\n"
687 " BRIK/HEAD (byte) files in a directory called PREFIX; \n"
688 " using AFNI_MAP is like using AFNI, but it gives non-\n"
689 " binarized *maps* of ROI connections.\n"
690 " Using BOTH produces AFNI and DUMP formats of outputs.\n"
691 " -dump_no_labtab :if the ROIS file has a label table, the default is to\n"
692 " use it in naming a '-dump_rois' output (if being used);\n"
693 " using this switch turn that off-- output file names \n"
694 " will be the same as if no label table were present.\n"
695 " -dump_lab_consec :if using `-dump_rois', then DON'T apply the numerical\n"
696 " labels of the original ROIs input to the output names.\n"
697 " This would only matter if input ROI labels aren't \n"
698 " consecutive and starting with one (e.g., if instead \n"
699 " they were 1,2,3,5,8,..).\n"
700 " ---> This is the opposite from previous default behavior, where\n"
701 " the option '-lab_orig_rois' was used to switch away\n"
702 " from consecutive-izing the labels in the output.\n"
703 " -posteriori :switch to have a bunch of individual files output, with\n"
704 " the value in each being the number of tracks per voxel\n"
705 " for that pair; works with '-dump_rois {AFNI | BOTH }',\n"
706 " where you get track-path maps instead of masks; makes\n"
707 " threshold for number of tracks between ROIs to keep to\n"
708 " be one automatically, regardless of setting in algopt.\n"
709 " -rec_orig :record dataset origin in the header of the *.trk file.\n"
710 " As of Sept. 2012, TrackVis doesn't use this info so it\n"
711 " wasn't included, but if you might want to map your \n"
712 " *.trk file later, then use the switch as the datasets's\n"
713 " Origin is necessary info for the mapping (the default\n"
714 " image in TrackVis might not pop up in the center of the\n"
715 " viewing window, though, just be aware). NB: including\n"
716 " the origin might become default at some point in time.\n"
717 " -do_trk_out :Switch ON outputting *.trk files, which are mainly to\n"
718 " be viewed in TrackVis (Wang et al., 2007). \n"
719 " (Feb, 2015): Default is to NOT output *.trk files.\n"
720 " -trk_opp_orient :If outputting *.trk files, you can choose to oppositize\n"
721 " the voxel_order parameter for the TRK file (only).\n"
722 " Thus, when inputting AFNI files with orient RAI, the\n"
723 " *.trk file would have voxel_order LPS; this is so that\n"
724 " files can be viewed in some other software, such as \n"
725 " DTK.\n"
726 " -nifti :output the PAIRMAP, INDIMAP, and any '-dump_rois' in\n"
727 " *.nii.gz format (default is BRIK/HEAD).\n"
728 " -no_indipair_out :Switch off outputting *INDIMAP* and *PAIRMAP* volumes.\n"
729 " This is probably just if you want to save file space;\n"
730 " also, for connectome-y studies with many (>100) target\n"
731 " regions, the output INDI and PAIR maps can be quite\n"
732 " large and/or difficult to write out. In some cases, it\n"
733 " might be better to just use '-dump_rois AFNI' instead.\n"
734 " Default is to output the INDI and PAIR map files.\n"
735 " -write_rois :write out a file (PREFIX.roi.labs) of all the ROI \n"
736 " (re-)labels, for example if the input ROIs aren't\n"
737 " simply consecutive and starting from 1. File has 3cols:\n"
738 " Input_ROI Condensed_form_ROI Power_of_2_label\n"
739 " -write_opts :write out all the option values into PREFIX.niml.opts.\n"
740 " -pair_out_power :switch to affect output of *PAIRMAP* output files. \n"
741 " Now, the default format is to output the >0 bricks with\n"
742 " tracks labelled by the target integers themselves.\n"
743 " This is not a unique labelling system, but it *is* far\n"
744 " easier to view and understand what's going on than\n"
745 " using a purely unique system based on using powers of\n"
746 " two of the ROIs (with integer summation for overlaps).\n"
747 " Using the switch '-pair_out_power' will change the\n"
748 " output of bricks [1] and higher to contain\n"
749 " information on connections stored as powers of two, so\n"
750 " that there is a unique decomposition in terms of\n"
751 " overlapped connections. However, it's *far* easier to\n"
752 " use '-dump_rois {AFNI | BOTH }' to get individual mask\n"
753 " files of the ROIs clearly (as well as annoying to need\n"
754 " to calculate powers of two simply to visualize the\n"
755 " connections.\n"
756 " -----> when considering this option, see the 'LABELTABLE'\n"
757 " description up above for how the labels work, with\n"
758 " or without an explicit table being entered.\n"
759 " -verb VERB :verbosity level, default is 0.\n"
760 "\n"
761 "****************************************************************************\n"
762 "\n"
763 "+ ALGOPT FILE EXAMPLES (note that different MODES have some different opts):\n"
764 " For '-mode {DET | MINIP}, the nicely readable NIML format of algopt file\n"
765 " would have a file name ending '.niml.opts' and contain something like the:\n"
766 " following seven lines:\n"
767 " <TRACK_opts\n"
768 " Thresh_FA=\"0.2\"\n"
769 " Thresh_ANG=\"60.000000\"\n"
770 " Thresh_Len=\"20.000000\"\n"
771 " Nseed_X=\"2\"\n"
772 " Nseed_Y=\"2\"\n"
773 " Nseed_Z=\"2\" />\n"
774 " The values above are actually all default values, and such a file would be\n"
775 " output using the '-write_opts' flag. For the same modes, one could get\n"
776 " the same result using a plain column of numbers, whose meaning is defined\n"
777 " by their order, contained in a file NOT ending in .niml.opts, such as \n"
778 " exemplified in the next six lines:\n"
779 " 0.2\n"
780 " 60\n"
781 " 20\n"
782 " 2\n"
783 " 2\n"
784 " 2\n"
785 " For '-mode PROB', the nice NIML format algopt file would contain something\n"
786 " like the next seven lines (again requiring the file name to end in\n"
787 " '.niml.opts'):\n"
788 " <TRACK_opts\n"
789 " Thresh_FA=\"0.2\"\n"
790 " Thresh_ANG=\"60.0\"\n"
791 " Thresh_Len=\"20.0\"\n"
792 " Thresh_Frac=\"0.001\"\n"
793 " Nseed_Vox=\"5\"\n"
794 " Nmonte=\"1000\" />\n"
795 " Again, those represent the default values, and could be given as a plain\n"
796 " column of numbers, in that order.\n"
797 "\n"
798 "* * ** * ** * ** * ** * ** * ** * ** * ** * ** * ** * ** * ** * ** * ** * **\n"
799 "\n"
800 "+ DTI LIST FILE EXAMPLE:\n"
801 " Consider, for example, if you hadn't used the '-sep_dsets' option when\n"
802 " outputting all the tensor information from 3dDWItoDT. Then one could\n"
803 " specify the DTI inputs for this program with a file called, e.g., \n"
804 " FILE_DTI_IN.niml.opts (the name *must* end with '.niml.opts'):\n"
805 " <DTIFILE_opts\n"
806 " dti_V1=\"SINGLEDT+orig[9..11]\"\n"
807 " dti_V2=\"SINGLEDT+orig[12..14]\"\n"
808 " dti_V3=\"SINGLEDT+orig[15..17]\"\n"
809 " dti_FA=\"SINGLEDT+orig[18]\"\n"
810 " dti_MD=\"SINGLEDT+orig[19]\"\n"
811 " dti_L1=\"SINGLEDT+orig[6]\"\n"
812 " dti_RD=\"SINGLEDT+orig[20]\" />\n"
813 " This represents the *minimum* set of input files needed when running\n"
814 " 3dTrackID. (Oct. 2016: RD now output by 3dDWItoDT, and not calc'ed \n"
815 " internally by 3dTrackID.)\n"
816 " One could also input extra data: an 'extra file' (XF) to take the place\n"
817 " of an FA map for determining where tracks can propagate; and up to four\n"
818 " other data sets (P1, P2, P3 and P4, standing for 'plus one' etc.) for\n"
819 " calculating mean/stdev properties in the obtained WM-ROIs:\n"
820 " <DTIFILE_opts\n"
821 " dti_V1=\"SINGLEDT+orig[9..11]\"\n"
822 " dti_V2=\"SINGLEDT+orig[12..14]\"\n"
823 " dti_V3=\"SINGLEDT+orig[15..17]\"\n"
824 " dti_XF=\"Segmented_WM.nii.gz\"\n"
825 " dti_FA=\"SINGLEDT+orig[18]\"\n"
826 " dti_MD=\"SINGLEDT+orig[19]\"\n"
827 " dti_L1=\"SINGLEDT+orig[6]\"\n"
828 " dti_RD=\"SINGLEDT+orig[20]\"\n"
829 " dti_P1=\"SINGLEDT+orig[7]\"\n"
830 " dti_P2=\"SINGLEDT+orig[8]\"\n"
831 " dti_P3=\"T1_map.nii.gz\"\n"
832 " dti_P4=\"PD_map.nii.gz\" />\n"
833 "\n"
834 "****************************************************************************\n"
835 "\n"
836 "+ EXAMPLES:\n"
837 " Here are just a few scenarios-- please see the Demo data set for *maaany*\n"
838 " more, as well as for fuller descriptions. To obtain the Demo, type the\n"
839 " following into a commandline:\n"
840 " $ @Install_FATCAT_demo\n"
841 " This will also unzip the archive, which contains required data (so it's\n"
842 " pretty big, currently >200MB), a README.txt file, and several premade\n"
843 " scripts that are ~heavily commented.\n"
844 "\n"
845 " A) Deterministic whole-brain tracking; set of targets is just the volume\n"
846 " mask. This can be useful for diagnostic purposes, sanity check for\n"
847 " gradients+data, for interactively selecting interesting subsets later,\n"
848 " etc. This uses most of the default algopts, but sets a higher minimum\n"
849 " length for keeping tracks:\n"
850 " $ 3dTrackID -mode DET \\\n"
851 " -dti_in DTI/DT \\\n"
852 " -netrois mask_DWI+orig \\\n"
853 " -logic OR \\\n"
854 " -alg_Thresh_Len 30 \\\n"
855 " -prefix DTI/o.WB\n"
856 "\n"
857 " B) Mini-probabilistic tracking through a multi-brik network file using a\n"
858 " DTI model and AND-logic. Instead of using the thresholded FA map to\n"
859 " guide tracking, an extra data set (e.g., a mapped anatomical\n"
860 " segmentation image) is input as the WM proxy; as such, what used to be\n"
861 " a threshold for adult parenchyma FA is now changed to an appropriate\n"
862 " value for the segmentation percentages; and this would most likely\n"
863 " also assume that 3dDWUncert had been used to calculate tensor value\n"
864 " uncertainties:\n"
865 " $ 3dTrackID -mode MINIP \\\n"
866 " -dti_in DTI/DT \\\n"
867 " -dti_extra T1_WM_in_DWI.nii.gz \\\n"
868 " -netrois ROI_ICMAP_GMI+orig \\\n"
869 " -logic AND \\\n"
870 " -mini_num 7 \\\n"
871 " -uncert DTI/o.UNCERT_UNC+orig. \\\n"
872 " -alg_Thresh_FA 0.95 \\\n"
873 " -prefix DTI/o.MP_AND_WM \n"
874 "\n"
875 " C) Full probabilistic tracking through a multi-brik network file using\n"
876 " HARDI-Qball reconstruction. The designated GFA file is used to guide\n"
877 " the tracking, with an appropriate threshold set and a smaller minimum\n"
878 " uncertainty of that GFA value (from this and example B, note how\n"
879 " generically the '-alg_Thresh_FA' functions, always setting a value for\n"
880 " the WM proxy map, whether it be literally FA, GFA or the dti_extra\n"
881 " file). Since HARDI-value uncertainty isn't yet calculable in AFNI,\n"
882 " brain-wide uniform values were assigned to the GFA and directions:\n"
883 " $ 3dTrackID -mode PROB \\\n"
884 " -hardi_gfa QBALL/GFA.nii.gz \\\n"
885 " -hardi_dirs QBALL/dirs.nii.gz \\\n"
886 " -netrois ROI_ICMAP_GMI+orig \\\n"
887 " -uncert QBALL/UNIFORM_UNC+orig. \\\n"
888 " -mask mask_DWI+orig \\\n"
889 " -alg_Thresh_FA 0.04 \\\n"
890 " -unc_min_FA 0.003 \\\n"
891 " -prefix QBALL/o.PR_QB \n"
892 "\n"
893 "****************************************************************************\n"
894 "\n"
895 " If you use this program, please reference the workhorse FACTID\n"
896 " tractography algorithm:\n"
897 " Taylor PA, Cho K-H, Lin C-P, Biswal BB (2012). Improving DTI\n"
898 " Tractography by including Diagonal Tract Propagation. PLoS ONE\n"
899 " 7(9): e43415. \n"
900 " and the introductory/description paper for FATCAT:\n"
901 " Taylor PA, Saad ZS (2013). FATCAT: (An Efficient) Functional And\n"
902 " Tractographic Connectivity Analysis Toolbox. Brain Connectivity.\n\n");
903 return;
904 }
905
906
907
main(int argc,char * argv[])908 int main(int argc, char *argv[])
909 {
910
911 TRACK_RUN_PARAMS InOpts;
912 int i;
913 int iarg;
914
915 // fill option struct with defaults
916 InOpts = SetInitialTrackingDefaults();
917
918 // ****************************************************************
919 // ****************************************************************
920 // load AFNI stuff
921 // ****************************************************************
922 // ****************************************************************
923
924 // INFO_message("version: CHI");
925 mainENTRY("3dTrackID"); machdep();
926
927 // scan args
928 if (argc == 1) { usage_TrackID(1); exit(0); }
929 iarg = 1;
930 while( iarg < argc && argv[iarg][0] == '-' ){
931 if( strcmp(argv[iarg],"-help") == 0 ||
932 strcmp(argv[iarg],"-h") == 0 ) {
933 usage_TrackID(strlen(argv[iarg])>3 ? 2:1);
934 exit(0);
935 }
936
937 if( strcmp(argv[iarg],"-verb") == 0) {
938 if( ++iarg >= argc )
939 ERROR_exit("Need argument after '-verb'") ;
940 set_tract_verb(atoi(argv[iarg]));
941 iarg++ ; continue ;
942 }
943
944 if( strcmp(argv[iarg],"-write_opts") == 0) {
945 InOpts.dump_opts=1;
946 iarg++ ; continue ;
947 }
948
949 if( strcmp(argv[iarg],"-write_rois") == 0) {
950 InOpts.ROIS_OUT=1;
951 iarg++ ; continue ;
952 }
953
954 if( strcmp(argv[iarg],"-pair_out_power") == 0) {
955 InOpts.PAIRPOWER=1;
956 iarg++ ; continue ;
957 }
958
959 if( strcmp(argv[iarg],"-pair_out_int") == 0) {
960 WARNING_message("(Harmless note:) the option '-pair_out_int' will"
961 " not be used anymore,\n"
962 "\t as it describes default behavior of PAIRMAP"
963 " file output;\n"
964 "\t see '-pair_out_power' for the ugly but unique"
965 " file labelling output.");
966 iarg++ ; continue ;
967 }
968
969 if( strcmp(argv[iarg],"-extra_tr_par") == 0) {
970 // scale NT grids by:
971 // 1) target vol,
972 // 2) target SA,
973 // 3) and target SA near WM.
974 InOpts.EXTRA_TR_PAR = N_XTR_PAR_LABS;
975 iarg++ ; continue ;
976 }
977
978 if( strcmp(argv[iarg],"-dump_lab_consec") == 0) {
979 InOpts.DUMP_ORIG_LABS=0;
980 iarg++ ; continue ;
981 }
982
983 if( strcmp(argv[iarg],"-nifti") == 0) {
984 InOpts.NIFTI_OUT=1;
985 iarg++ ; continue ;
986 }
987
988 if( strcmp(argv[iarg],"-netrois") == 0 ){
989 if( ++iarg >= argc )
990 ERROR_exit("Need argument after '-netrois'") ;
991
992 InOpts.NAMEIN_netrois = strdup(argv[iarg]);
993
994 iarg++ ; continue ;
995 }
996
997 if( strcmp(argv[iarg],"-uncert") == 0 ){
998 if( ++iarg >= argc )
999 ERROR_exit("Need argument after '-uncert'") ;
1000
1001 InOpts.NAMEIN_uncert = strdup(argv[iarg]);
1002
1003 iarg++ ; continue ;
1004 }
1005
1006 if( strcmp(argv[iarg],"-unc_min_FA") == 0) {
1007 if( ++iarg >= argc )
1008 ERROR_exit("Need argument after '-unc_min_FA'") ;
1009
1010 InOpts.NAMEIN_uncFA = strdup(argv[iarg]);
1011
1012 iarg++ ; continue ;
1013 }
1014
1015 if( strcmp(argv[iarg],"-unc_min_V") == 0) {
1016 if( ++iarg >= argc )
1017 ERROR_exit("Need argument after '-unc_min_V'") ;
1018
1019 InOpts.NAMEIN_uncEI = strdup(argv[iarg]);
1020
1021 iarg++ ; continue ;
1022 }
1023
1024 if( strcmp(argv[iarg],"-prefix") == 0 ){
1025 iarg++ ; if( iarg >= argc )
1026 ERROR_exit("Need argument after '-prefix'");
1027 InOpts.prefix = strdup(argv[iarg]) ;
1028 if( !THD_filename_ok(InOpts.prefix) )
1029 ERROR_exit("Illegal name after '-prefix'");
1030 iarg++ ; continue ;
1031 }
1032
1033 if( strcmp(argv[iarg],"-dti_in") == 0 ){
1034 iarg++ ; if( iarg >= argc )
1035 ERROR_exit("Need argument after '-in_dti'");
1036 InOpts.dti_inpref = strdup(argv[iarg]) ;
1037
1038 iarg++ ; continue ;
1039 }
1040
1041 if( strcmp(argv[iarg],"-dti_search_NO") == 0) {
1042 InOpts.NO_NONDTI_SEARCH = 1;
1043 iarg++ ; continue ;
1044 }
1045
1046 // tracts are chopped at targ surface bounds
1047 if( strcmp(argv[iarg],"-targ_surf_stop") == 0) {
1048 InOpts.TARG_SURF_STOP = 1;
1049 iarg++ ; continue ;
1050 }
1051 // ... or prevent from penetrating bounds
1052 if( strcmp(argv[iarg],"-targ_surf_twixt") == 0) {
1053 InOpts.TARG_SURF_STOP = 1;
1054 InOpts.TARG_SURF_TWIXT = 1;
1055 iarg++ ; continue ;
1056 }
1057
1058 if( strcmp(argv[iarg],"-dti_list") == 0 ){
1059 iarg++ ;
1060 if( iarg >= argc )
1061 ERROR_exit("Need argument after '-dti_list'");
1062 InOpts.dti_listname = strdup(argv[iarg]) ;
1063
1064 iarg++ ; continue ;
1065 }
1066
1067 if( strcmp(argv[iarg],"-algopt") == 0 ){
1068 iarg++ ;
1069 if( iarg >= argc )
1070 ERROR_exit("Need argument after '-algopt'");
1071 InOpts.algopt_file_name = strdup(argv[iarg]) ;
1072
1073 iarg++ ; continue ;
1074 }
1075
1076 // can also set options at commandline
1077 if( strcmp(argv[iarg],"-alg_Thresh_FA") == 0) {
1078 iarg++ ;
1079 if( iarg >= argc )
1080 ERROR_exit("Need argument after '-alg_Thresh_FA'");
1081
1082 InOpts.MinFA = atof(argv[iarg]);
1083
1084 iarg++ ; continue ;
1085 }
1086 if( strcmp(argv[iarg],"-alg_Thresh_ANG") == 0) {
1087 iarg++ ;
1088 if( iarg >= argc )
1089 ERROR_exit("Need argument after '-alg_Thresh_Ang'");
1090 InOpts.MaxAngDeg = atof(argv[iarg]);
1091
1092 iarg++ ; continue ;
1093 }
1094 if( strcmp(argv[iarg],"-alg_Thresh_Len") == 0) {
1095 iarg++ ;
1096 if( iarg >= argc )
1097 ERROR_exit("Need argument after '-alg_Thresh_Len'");
1098 InOpts.MinL = atof(argv[iarg]);
1099 iarg++ ; continue ;
1100 }
1101 if( strcmp(argv[iarg],"-alg_Thresh_Frac") == 0) {
1102 iarg++ ;
1103 if( iarg >= argc )
1104 ERROR_exit("Need argument after '-alg_Thresh_Frac'");
1105 InOpts.NmNsFr = atof(argv[iarg]);
1106
1107 iarg++ ; continue ;
1108 }
1109 if( strcmp(argv[iarg],"-alg_Nseed_Vox") == 0) {
1110 iarg++ ;
1111 if( iarg >= argc )
1112 ERROR_exit("Need argument after '-alg_Nseed_Vox'");
1113 InOpts.Nseed = atoi(argv[iarg]);
1114
1115 iarg++ ; continue ;
1116 }
1117 if( strcmp(argv[iarg],"-alg_Nmonte") == 0) {
1118 iarg++ ;
1119 if( iarg >= argc )
1120 ERROR_exit("Need argument after '-alg_Nmonte'");
1121 InOpts.Nmonte = atoi(argv[iarg]);
1122
1123 iarg++ ; continue ;
1124 }
1125 if( strcmp(argv[iarg],"-alg_Nseed_X") == 0) {
1126 iarg++ ;
1127 if( iarg >= argc )
1128 ERROR_exit("Need argument after '-alg_Nseed_X'");
1129 InOpts.SeedPerV[0] = atoi(argv[iarg]);
1130
1131 iarg++ ; continue ;
1132 }
1133 if( strcmp(argv[iarg],"-alg_Nseed_Y") == 0) {
1134 iarg++ ;
1135 if( iarg >= argc )
1136 ERROR_exit("Need argument after '-alg_Nseed_Y'");
1137 InOpts.SeedPerV[1] = atoi(argv[iarg]);
1138
1139 iarg++ ; continue ;
1140 }
1141 if( strcmp(argv[iarg],"-alg_Nseed_Z") == 0) {
1142 iarg++ ;
1143 if( iarg >= argc )
1144 ERROR_exit("Need argument after '-alg_Nseed_Z'");
1145 InOpts.SeedPerV[2] = atoi(argv[iarg]);
1146
1147 iarg++ ; continue ;
1148 }
1149
1150 // Not in main code for now
1151 if( strcmp(argv[iarg],"-choose_seed") == 0) {
1152 iarg++ ;
1153 if( iarg >= argc )
1154 ERROR_exit("Need argument after '-choose_seed'");
1155 InOpts.CHOOSE_SEED = atoi(argv[iarg]);
1156
1157 iarg++ ; continue ;
1158 }
1159
1160 if( strcmp(argv[iarg],"-no_trk_out") == 0) {
1161 WARNING_message("The switch '-no_trk_out' is no longer used."
1162 " As of Feb, 2015, TrackVis-style *.trk"
1163 "\n files are no longer created by default."
1164 " To output them, use '-do_trk_out'." );
1165 iarg++ ; continue ;
1166 }
1167 if( strcmp(argv[iarg],"-do_trk_out") == 0) {
1168 InOpts.OUTPUT_TRK = 1;
1169 iarg++ ; continue ;
1170 }
1171 if( strcmp(argv[iarg],"-trk_opp_orient") == 0) {
1172 InOpts.TRK_OPP_ORI = 1;
1173 iarg++ ; continue ;
1174 }
1175
1176 if( strcmp(argv[iarg],"-no_indipair_out") == 0) {
1177 InOpts.OUT_INDIPAIR = 0;
1178 iarg++ ; continue ;
1179 }
1180
1181 if( strcmp(argv[iarg],"-dti_extra") == 0) {
1182 if( ++iarg >= argc )
1183 ERROR_exit("Need argument after '-dti_extra'");
1184
1185 InOpts.in_EXTRA = strdup(argv[iarg]) ;
1186
1187 iarg++ ; continue ;
1188 }
1189
1190 if( strcmp(argv[iarg],"-mask") == 0 ){
1191 iarg++ ; if( iarg >= argc )
1192 ERROR_exit("Need argument after '-mask'");
1193
1194 InOpts.NAMEIN_mask = strdup( argv[iarg]);
1195
1196 iarg++ ; continue ;
1197 }
1198
1199 if( strcmp(argv[iarg],"-thru_mask") == 0 ){
1200 iarg++ ; if( iarg >= argc )
1201 ERROR_exit("Need argument after '-thru_mask'");
1202
1203 InOpts.NAMEIN_thru = strdup( argv[iarg]);
1204
1205 iarg++ ; continue ;
1206 }
1207
1208
1209 if( strcmp(argv[iarg],"-uncut_at_rois") == 0) {
1210 InOpts.ONLY_BT=0;
1211 iarg++ ; continue ;
1212 }
1213
1214 if( strcmp(argv[iarg],"-posteriori") == 0) {
1215 InOpts.POST=1;
1216 iarg++ ; continue ;
1217 }
1218
1219 if( strcmp(argv[iarg],"-mode") == 0) { // @)
1220 iarg++ ;
1221 if( iarg >= argc )
1222 ERROR_exit("Need argument after '-mode'");
1223
1224 InOpts.NAMEIN_mode = strdup(argv[iarg]);
1225
1226 iarg++ ; continue ;
1227 }
1228
1229 if( strcmp(argv[iarg],"-logic") == 0) { // @)
1230 iarg++ ; if( iarg >= argc )
1231 ERROR_exit("Need argument after '-logic'");
1232
1233 INFO_message("ROI logic type is: %s",argv[iarg]);
1234 if( strcmp(argv[iarg],"AND") == 0 )
1235 InOpts.LOG_TYPE = 1;
1236 else if( strcmp(argv[iarg],"OR") == 0 )
1237 InOpts.LOG_TYPE = 0;
1238 else
1239 ERROR_exit("Illegal after '-logic': need 'OR' or 'AND'");
1240
1241 iarg++ ; continue ;
1242 }
1243
1244 if( strcmp(argv[iarg],"-hardi_dirs") == 0) { // @)
1245 iarg++ ; if( iarg >= argc )
1246 ERROR_exit("Need argument after '-hardi_dirs'");
1247 InOpts.hardi_dir = strdup(argv[iarg]) ;
1248
1249 iarg++ ; continue ;
1250 }
1251
1252 if( strcmp(argv[iarg],"-hardi_gfa") == 0) { // @)
1253 iarg++ ; if( iarg >= argc )
1254 ERROR_exit("Need argument after '-hardi_gfa'");
1255 InOpts.hardi_gfa = strdup(argv[iarg]) ;
1256
1257 iarg++ ; continue ;
1258 }
1259
1260 if( strcmp(argv[iarg],"-hardi_pars") == 0) { // @)
1261 iarg++ ; if( iarg >= argc )
1262 ERROR_exit("Need argument after '-hardi_pars'");
1263 InOpts.hardi_pars = strdup(argv[iarg]) ;
1264
1265 iarg++ ; continue ;
1266 }
1267
1268 // can do a mini monte carlo if you want; for example, use 1
1269 // seed per ROI and then see more locs to tracks possibly; prob
1270 // don't want to make tooo large for visualization stuff.
1271 if( strcmp(argv[iarg],"-mini_num") == 0) { // @)
1272 iarg++ ; if( iarg >= argc )
1273 ERROR_exit("Need argument after '-mini_num'");
1274
1275 InOpts.MINI_PROB_NM = atoi(argv[iarg]);
1276
1277 iarg++ ; continue ;
1278 }
1279
1280 // Feb,2015
1281 if( strcmp(argv[iarg],"-bundle_thr") == 0) { // @)
1282 iarg++ ; if( iarg >= argc )
1283 ERROR_exit("Need argument after '-bundle_thr'");
1284
1285 InOpts.THRESH_BUNDS = atoi(argv[iarg]);
1286
1287 iarg++ ; continue ;
1288 }
1289
1290
1291 /* TURNED OFF FOR MOMENT, WORK IN PROGRESS
1292 if( strcmp(argv[iarg],"-jump_les") == 0) { // @)) prop check
1293 JUMPLES=1;
1294 iarg++ ; continue ;
1295 } */
1296
1297 if( strcmp(argv[iarg],"-rec_orig") == 0) { // @)
1298 InOpts.RECORD_ORIG=1;
1299 iarg++ ; continue ;
1300 }
1301
1302 if( strcmp(argv[iarg],"-tract_out_mode") == 0) { // @)
1303 if( ++iarg >= argc )
1304 ERROR_exit("Need argument after '-tract_out_mode'") ;
1305
1306 InOpts.NAMEIN_outmode = strdup(argv[iarg]);
1307
1308 iarg++ ; continue ;
1309 }
1310
1311 if( strcmp(argv[iarg],"-dump_rois") == 0) {
1312 iarg++ ; if( iarg >= argc )
1313 ERROR_exit("Need argument after '-dump_rois'");
1314
1315 if( strcmp(argv[iarg],"DUMP") == 0 )
1316 InOpts.DUMP_TYPE = 1;
1317 else if( strcmp(argv[iarg],"AFNI") == 0 )
1318 InOpts.DUMP_TYPE = 2;
1319 else if( strcmp(argv[iarg],"AFNI_MAP") == 0 )
1320 InOpts.DUMP_TYPE = 4;
1321 else if( strcmp(argv[iarg],"BOTH") == 0 )
1322 InOpts.DUMP_TYPE = 3;
1323 else
1324 ERROR_exit("Illegal after '-dump_rois': need 'DUMP',"
1325 " 'AFNI' or 'BOTH' or the new 'AFNI_MAP'.");
1326
1327 iarg++ ; continue ;
1328 }
1329
1330 if( strcmp(argv[iarg],"-dump_no_labtab") == 0) {
1331 // if NETROIS has a label table, default is to use it in
1332 // naming a dump_rois output; this can turn that off.
1333 InOpts.DUMP_with_LABELS=0;
1334 iarg++ ; continue ;
1335 }
1336
1337 if( strcmp(argv[iarg],"-cut_at_rois") == 0) {
1338 ERROR_exit("No longer used! Default behavior now to cut_at_rois;"
1339 " turn *off* with '-uncut_at_rois', if desired");
1340 iarg++ ; continue ;
1341 }
1342
1343 if( strcmp(argv[iarg],"-det_net") == 0) {
1344 ERROR_exit("No longer used! use '-logic {AND | OR}', instead.");
1345 iarg++ ; continue ;
1346 }
1347
1348 ERROR_message("Bad option '%s'\n",argv[iarg]) ;
1349 suggest_best_prog_option(argv[0], argv[iarg]);
1350 exit(1);
1351 }
1352
1353 if (iarg < 5) {
1354 ERROR_message("Too few options. Try -help for details.\n");
1355 exit(1);
1356 }
1357
1358 // DONE FILLING, now call
1359 i = RunTrackingMaestro(1, InOpts, argc, argv);
1360
1361 if( !i )
1362 ERROR_message("Error in run?");
1363
1364 exit(0);
1365 }
1366
1367
1368 // BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
1369
1370
RunTrackingMaestro(int comline,TRACK_RUN_PARAMS opts,int argc,char * argv[])1371 int RunTrackingMaestro( int comline, TRACK_RUN_PARAMS opts,
1372 int argc, char *argv[])
1373 {
1374
1375 int i,j,k,m,n,ii,jj,kk,mm,gg,nn,hh,bb,cc,rr,ss,uu,vv,ll,lll;
1376 int nvox_rois=0;
1377 int nvox_unc=0;
1378 THD_3dim_dataset *insetUC = NULL;
1379 THD_3dim_dataset **insetPARS=NULL; // DEF_DTI: [0] X, [1] FA, [2] MD, [3] L1
1380 THD_3dim_dataset **insetVECS=NULL; // DEF_DTI: [0] V1, [1] V2, [2] V3
1381 MRI_IMAGE *temp_vec=NULL;
1382
1383 // default for DTI, but easily changed for HARDI
1384 int PARS_BOT = 1; // default is to start at brik 1 for pars, with FA there.
1385 int PARS_TOP = 0; // number of bricks.
1386 int PARS_N = 0;
1387 int Noutmat = 0; // number of output matrices in GRID file
1388 char **ParLab=NULL;
1389 int DEF_DTI = 1; // default mode, mainly for # of params
1390 char **wild_names=NULL; //[N_DTI_MAX_PARS][32]; // get names
1391 int pref_offset = 0;
1392 char temp_name[32];
1393 int hardi_pref_len=0;
1394 int FOUND_DTI;
1395 int tot_FOUND_DTI=0;
1396 int TR_MODE = -1; // -1=ERR; 0=DET; 1=MINIP; 2=PROB;
1397 char *list_modes[3] = {"DET", "MINIP", "PROB"}; // match with PARLABS later
1398
1399 THD_3dim_dataset *ROI_set=NULL;
1400
1401 char ***gdset_roi_names=NULL;
1402 SUMA_DSET *gset=NULL;
1403 float ***flat_matr=NULL;
1404 float *xyz=NULL;
1405 char OUT_gdset[THD_MAX_NAME];
1406 NI_group *GDSET_netngrlink=NULL;
1407 char *NAME_gdset=NULL;
1408
1409 THD_3dim_dataset *thrumask=NULL;
1410 int doesnt_pass;
1411 int KEEP_GOING=1;
1412
1413 THD_3dim_dataset *MASK=NULL;
1414 int HAVE_MASK=0;
1415 int ***mskd=NULL; // define mask of where time series are nonzero
1416 int ni=0;
1417 time_t t_start;
1418
1419 //char *DTI_SCAL_INS[N_dti_scal] = {"FA", "MD", "L1"}; // match PARLABS later
1420 //char *DTI_VECT_INS[N_dti_vect] = {"V1", "V2", "V3"};
1421 //char *DTI_PLUS_INS[N_plus_dtifile] = {"P1", "P2", "P3", "P4"};
1422
1423 int EXTRAFILE=0; // switch for whether other file is input as WM map
1424
1425 char OUT_grid[THD_MAX_NAME];
1426 char OUT_map[THD_MAX_NAME];
1427 char OUT_mask[THD_MAX_NAME];
1428 char OUT_rois[THD_MAX_NAME];
1429 FILE *fin4, *fout1;
1430
1431 float MaxAng=0.;
1432 int NmNsThr=0;
1433 int ArrMax=0;
1434 float unc_minfa_std=0.015, unc_minei_std=0.06; // min uncert stds.
1435 int USER_UNC_FA=0, USER_UNC_EI=0; // whether user has input explicit vals
1436
1437 int Nvox=-1; // tot number vox
1438 int Ndata=-1;
1439 float VoxVol=0;
1440 int *Dim=NULL; //[3]={0,0,0}; // dim in each dir
1441 float **LocSeed=NULL; // fractional locations within voxel,
1442 float Ledge[3]; // voxel edge lengths
1443
1444 int **Tforw=NULL, **Tback=NULL;
1445 int **Ttot=NULL;
1446 float **flTforw=NULL, **flTback=NULL;
1447 float **coorded=NULL;
1448 float **copy_coorded=NULL; // copy will have perturbed info;
1449 float **UNC=NULL; // for DTI now: just combined pieces, not bias+stdev
1450 int ***INDEX=NULL,***INDEX2=NULL;
1451 int **MAPROI=NULL; // store what ROIs are where, per data voxel
1452 int ***NETROI=NULL; // store connection info
1453 int len_forw, len_back; // int count of num of squares through
1454 float phys_forw[1], phys_back[1];
1455 int idx,idx3;
1456
1457 int in[3]; // to pass to trackit
1458 float physin[3]; // also for trackit, physical loc,
1459 int totlen;
1460 float totlen_phys, phys_l;
1461 int Numtract;
1462
1463 int READS_in;
1464 float READS_fl;
1465 int onoff[3];
1466 int BreakAddCont=0;
1467
1468 char dset_or[4] = "RAI";
1469 char *voxel_order=NULL;//[4]="---";
1470 THD_3dim_dataset *dsetn=NULL;
1471 int *TV_switch=NULL;//[3] = {0,0,0};
1472
1473 int **Prob_grid=NULL; // NROIx(NROI+1)/2 grid (NROI is num of ROIs);
1474 float ***Prob_grid_L=NULL; // for bundle lengths; Nov,2016
1475 float ***Param_grid=NULL; // same size as Prob_Grid (extra dim for
1476 // ROI stats)
1477 int *list_rois=NULL, *temp_list=NULL;
1478 int MAXNROI=0;
1479 int N_nets=0;
1480
1481 float w2,w3,thetval;
1482 float tempv[3];
1483 float tempvmagn;
1484 float testang;
1485
1486 int *NROI=NULL,*INVROI=NULL;
1487 int **ROI_LABELS=NULL, **INV_LABELS=NULL; // allow labels to be non-consec
1488 NI_element *nel=NULL;
1489 char *postfix[4]={"+orig.HEAD\0",".nii.gz\0",".nii\0","+tlrc.HEAD\0"};
1490 int FOUND = -1;
1491
1492 // @) this block of values for detnet options
1493 int DETNET = 0; // switch on detnet
1494 TAYLOR_NETWORK **tnet=NULL;
1495 TAYLOR_BUNDLE ***tb=NULL; // @))), for new IO
1496 int *N_bund=NULL; // @))) number of possible tri+diag bund/net
1497 TAYLOR_TRACT **tt=NULL;
1498 char *mode = "NI_fast_binary";
1499 int *id=NULL;
1500 float **flTtot=NULL;
1501 int start_loc=0;
1502 int trL=0;
1503 char **prefix_det=NULL;
1504 char OUT_bin[THD_MAX_NAME];
1505 FILE **fout0=NULL;
1506 tv_io_header headerDTI = {.id_string = "TRACK\0",
1507 .origin = {0,0,0},
1508 .n_scalars = 3, // eventually can adjust;
1509 .scal_n[0] = "FA",
1510 .scal_n[1] = "MD",
1511 .scal_n[2] = "L1",
1512 .n_properties = 0,
1513 .vox_to_ras = {{0.,0.,0.,0.},{0.,0.,0.,0.},
1514 {0.,0.,0.,0.},{0.,0.,0.,0.}},
1515 // reset this later based on actual data set
1516 .voxel_order = "RAI\0",
1517 .invert_x = 0,
1518 .invert_y = 0,
1519 .invert_z = 0,
1520 .swap_xy = 0,
1521 .swap_yz = 0,
1522 .swap_zx = 0,
1523 .n_count = 0,
1524 .version = 2,
1525 .hdr_size = 1000};
1526 tv_io_header headerHAR = {.id_string = "TRACK\0",
1527 .origin = {0,0,0},
1528 .n_scalars = 1, // niml has all data anyways...
1529 //.scal_n[0] = "FA",
1530 //.scal_n[1] = "MD",
1531 //.scal_n[2] = "L1",
1532 .n_properties = 0,
1533 .vox_to_ras = {{0.,0.,0.,0.},{0.,0.,0.,0.},
1534 {0.,0.,0.,0.},{0.,0.,0.,0.}},
1535 // reset this later based on actual data set
1536 .voxel_order = "RAI\0",
1537 .invert_x = 0,
1538 .invert_y = 0,
1539 .invert_z = 0,
1540 .swap_xy = 0,
1541 .swap_yz = 0,
1542 .swap_zx = 0,
1543 .n_count = 0,
1544 .version = 2,
1545 .hdr_size = 1000};
1546
1547 float Orig[3] = {0.0,0.0,0.0};
1548
1549 int vA,vA0,vB,vB0,vB1; // will be indices along tracks in NOT-ROI checking
1550 int TESTCOUNT=0;
1551
1552 int N_HAR = 0;
1553 int N_uncert = 0;
1554 short *DirPerVox=NULL; // can have diff num of dir/vox
1555 int N_HARPARS=0;
1556
1557 // @)) jump lesion searching stuff, to be used prob only with DETNET
1558 int JUMPLES=0; // switch to turn on prop check.
1559 float ***TROUT=NULL;
1560
1561 char *wild_list=NULL;
1562 char **wglob=NULL, **wsort=NULL;
1563 int nglob, nsort, *isrt=NULL;
1564 int wild_all_files = 0, wild_ci=0;
1565 int wild_noext=2; // "-wild_files_noAext_noAview"
1566 int wild_orig_name = 1; // "-wild_files_orig_name"
1567
1568 char tprefix[THD_MAX_PREFIX];
1569 int nfiles;
1570
1571 int ***ROI_SIZES=NULL;
1572
1573 Dtable *roi_dtable=NULL;
1574 char *LabTabStr=NULL;
1575 char ***ROI_STR_LABELS=NULL;
1576 char EleNameStr[128];
1577
1578 // for random number generation
1579 const gsl_rng_type * T=NULL;
1580 gsl_rng *r=NULL;
1581 long seed1, seed2;
1582
1583
1584 // * * * * * * * * * * * * * * Load general stuff in * * * * * * * * * *
1585 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
1586
1587 // Use random seeds, unless special repeat option is on
1588 if( opts.CHOOSE_SEED )
1589 seed1 = seed2 = opts.CHOOSE_SEED;
1590 else {
1591 seed1 = time(NULL);
1592 seed2 = time(NULL) ;
1593 }
1594
1595 srand(seed1);
1596 gsl_rng_env_setup();
1597 T = gsl_rng_default;
1598 r = gsl_rng_alloc (T);
1599 gsl_rng_set (r, seed2);
1600
1601 // from '-mode'
1602 if( opts.NAMEIN_mode ) {
1603 INFO_message("Tracking mode: %s",opts.NAMEIN_mode);
1604 if( strcmp(opts.NAMEIN_mode,list_modes[0]) == 0 ) { // det
1605 DETNET = 1;
1606 TR_MODE = 0;
1607 }
1608 else if( strcmp(opts.NAMEIN_mode,list_modes[1]) == 0 ) { // minip
1609 DETNET = 1;
1610 TR_MODE = 1;
1611 }
1612 else if( strcmp(opts.NAMEIN_mode,list_modes[2]) == 0 ) { //prob
1613 DETNET = 0;
1614 TR_MODE = 2;
1615 }
1616 else
1617 ERROR_exit("Illegal after '-mode': need '%s' or '%s' or '%s'",
1618 list_modes[0], list_modes[1], list_modes[2]);
1619 }
1620 else {
1621 ERROR_message("Failed to load '-mode' option correctly.");
1622 exit(3);
1623 }
1624
1625 if( opts.hardi_dir || opts.hardi_gfa || opts.hardi_pars )
1626 DEF_DTI = 0;
1627
1628 if( opts.hardi_gfa )
1629 PARS_BOT = 0;
1630
1631 // from '-netrois'
1632 if( opts.NAMEIN_netrois ) {
1633 ROI_set = THD_open_dataset( opts.NAMEIN_netrois );
1634 if( ROI_set == NULL )
1635 ERROR_exit("Can't open netrois dataset '%s'", opts.NAMEIN_netrois);
1636 DSET_load(ROI_set) ; CHECK_LOAD_ERROR(ROI_set) ;
1637 nvox_rois = DSET_NVOX(ROI_set) ;
1638 N_nets = DSET_NVALS(ROI_set) ;
1639
1640 NROI = (int *)calloc(N_nets, sizeof(int));
1641 INVROI = (int *)calloc(N_nets, sizeof(int));
1642
1643 if( (NROI == NULL) || (INVROI == NULL) ) {
1644 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
1645 exit(122);
1646 }
1647
1648 for( i=0 ; i<N_nets ; i++)
1649 INVROI[i] = (int) THD_subbrick_max(ROI_set, i, 1);
1650
1651 ROI_LABELS = calloc( N_nets,sizeof(ROI_LABELS));
1652 for(i=0 ; i<N_nets ; i++)
1653 ROI_LABELS[i] = calloc(INVROI[i]+1,sizeof(int));
1654
1655 INV_LABELS = calloc( N_nets,sizeof(INV_LABELS));
1656 for(i=0 ; i<N_nets ; i++)
1657 INV_LABELS[i] = calloc(INVROI[i]+1,sizeof(int));
1658
1659 if( (ROI_LABELS == NULL) || (ROI_LABELS == NULL) ) {
1660 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
1661 exit(123);
1662 }
1663
1664 // after this, ROI_LABELS will have network labels store compactly
1665 // and NROI will have the actual number of ROI per netw, not just
1666 // the max label value
1667 // INV_LABELS has compact index labels at input label locations,
1668 // hence is inverse of ROI_LABELS
1669 // ** This function allows us to switch to talk about with i-th index
1670 // label is used, using NROI to define sizes of things
1671 bb = ViveLeRoi(ROI_set, ROI_LABELS, INV_LABELS, NROI, INVROI);
1672 if( bb != 1)
1673 ERROR_exit("Problem loading/assigning ROI labels");
1674
1675 for( i=0 ; i<N_nets ; i++) {
1676 if( NROI[i]>MAXNROI )
1677 MAXNROI = NROI[i];
1678 }
1679
1680 fprintf(stdout,"\n");
1681 for( bb=0 ; bb<N_nets ; bb++)
1682 INFO_message("Number of ROIs in netw[%d] = %d",bb,NROI[bb]);
1683
1684 ROI_STR_LABELS = (char ***) calloc( N_nets, sizeof(char **) );
1685 for ( i=0 ; i<N_nets ; i++ )
1686 ROI_STR_LABELS[i] = (char **) calloc( NROI[i]+1, sizeof(char *) );
1687 for ( i=0 ; i<N_nets ; i++ )
1688 for ( j=0 ; j<NROI[i]+1 ; j++ )
1689 ROI_STR_LABELS[i][j] = (char *) calloc( 100 , sizeof(char) );
1690 if( ROI_STR_LABELS == NULL ) {
1691 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
1692 exit(123);
1693 }
1694
1695 // Sept 2014: Labeltable stuff
1696 if ( ROI_set->Label_Dtable = DSET_Label_Dtable(ROI_set) ) {
1697 if ( LabTabStr = Dtable_to_nimlstring( DSET_Label_Dtable(ROI_set),
1698 "VALUE_LABEL_DTABLE") ) {
1699 //fprintf(stdout,"%s", LabTabStr);
1700 if (!(roi_dtable = Dtable_from_nimlstring(LabTabStr))) {
1701 ERROR_exit("Could not parse labeltable.");
1702 }
1703 }
1704 else {
1705 INFO_message("No label table from '-netrois'.");
1706 }
1707 }
1708
1709 bb = Make_ROI_Output_Labels( ROI_STR_LABELS,
1710 ROI_LABELS,
1711 N_nets,
1712 NROI,
1713 roi_dtable,
1714 opts.DUMP_with_LABELS);
1715 }
1716 else {
1717 ERROR_message("Failed to load '-netrois' options.");
1718 exit(2);
1719 }
1720
1721 // from '-uncert'
1722 if( opts.NAMEIN_uncert ) {
1723 insetUC = THD_open_dataset( opts.NAMEIN_uncert );
1724 if( insetUC == NULL )
1725 ERROR_exit("Can't open uncert dataset '%s'",
1726 opts.NAMEIN_uncert) ;
1727 DSET_load(insetUC) ; CHECK_LOAD_ERROR(insetUC);
1728 nvox_unc = DSET_NVOX(insetUC);
1729 N_uncert = DSET_NVALS(insetUC); // needed if HARDI input now
1730 }
1731
1732 if( opts.NAMEIN_uncFA ){
1733 unc_minfa_std = atof(opts.NAMEIN_uncFA);
1734 USER_UNC_FA = 1;
1735 }
1736
1737 if( opts.NAMEIN_uncEI ){
1738 unc_minei_std = atof(opts.NAMEIN_uncEI);
1739 unc_minei_std*= CONV;
1740 USER_UNC_EI = 1;
1741 }
1742
1743 // for naming stuff later
1744 if( opts.dti_inpref )
1745 hardi_pref_len = strlen(opts.dti_inpref); // assume an underscore
1746 else if( opts.hardi_pars )
1747 hardi_pref_len = strlen(opts.hardi_pars);
1748
1749 if( opts.NAMEIN_mask ){
1750 MASK = THD_open_dataset(opts.NAMEIN_mask) ;
1751 if( MASK == NULL )
1752 ERROR_exit("Can't open time series dataset '%s'.",opts.NAMEIN_mask);
1753
1754 DSET_load(MASK); CHECK_LOAD_ERROR(MASK);
1755 HAVE_MASK=1; // kind of needless switch...
1756 }
1757
1758 if( opts.in_EXTRA ){
1759 EXTRAFILE = 1; // switch on
1760 PARS_BOT = 0;
1761 }
1762
1763 if( opts.NAMEIN_outmode ) {
1764 if (strcmp(opts.NAMEIN_outmode, "NI_fast_binary") &&
1765 strcmp(opts.NAMEIN_outmode, "NI_fast_text") &&
1766 strcmp(opts.NAMEIN_outmode, "NI_slow_binary") &&
1767 strcmp(opts.NAMEIN_outmode, "NI_slow_text") ) {
1768 ERROR_message("Bad value (%s) for -tract_out_mode",
1769 opts.NAMEIN_outmode);
1770 exit(1);
1771 }
1772 mode = opts.NAMEIN_outmode;//argv[iarg];
1773 }
1774
1775 // * * * * * * * * * * * * * * check HARDI vs DTI * * * * * * * * * * * *
1776
1777 if( opts.algopt_file_name ) {
1778
1779 if( (TR_MODE == 0) || (TR_MODE == 1) ) { // DET or MINIP
1780 if (!(nel = ReadTractAlgOpts_M(opts.algopt_file_name))) {
1781 ERROR_message("Failed to read options in %s\n",
1782 opts.algopt_file_name);
1783 exit(19);
1784 }
1785
1786 if (NI_getTractAlgOpts_M(nel, &(opts.MinFA),&(opts.MaxAngDeg),
1787 &(opts.MinL),
1788 opts.SeedPerV)) {
1789 ERROR_message("Failed to get options");
1790 exit(1);
1791 }
1792 NI_free_element(nel); nel=NULL;
1793
1794 if( (opts.MinFA<0.001) &&
1795 (opts.MaxAngDeg<0.001) &&
1796 (opts.MinL<0.001) )
1797 ERROR_exit("All of alg_MinFA, alg_MaxAngDeg and alg_MinL"
1798 " are <0.001. This looks suspicious. Is your\n"
1799 "'-algopt' file %s written in the niml.opts"
1800 " format (i.e., with nice labels per parameter)\n"
1801 "but it doesn't end with .niml.opts, so it's being"
1802 " treated like the old school, plain-numbers-in-a-\n"
1803 "column format, and thus causing an error here?",
1804 opts.algopt_file_name);
1805
1806 if( (opts.SeedPerV[0] < 1) ||
1807 (opts.SeedPerV[1] < 1) ||
1808 (opts.SeedPerV[2] < 1) ) {
1809 WARNING_message("One of the dimensions of seeds was <1\n"
1810 "Perhaps you input a probabilistic algopt file?\n"
1811 "Since you are chose the deterministic switch, I\n"
1812 "am going to use 1 seed per vox.");
1813 opts.SeedPerV[0] = opts.SeedPerV[1] = opts.SeedPerV[2] = 1;
1814 }
1815
1816 }
1817 else { // PROB
1818 if (!(nel = ReadProbTractAlgOpts_M(opts.algopt_file_name))) {
1819 ERROR_message("Failed to read options in %s\n",
1820 opts.algopt_file_name);
1821 exit(19);
1822 }
1823
1824 if (NI_getProbTractAlgOpts_M(nel, &(opts.MinFA),&(opts.MaxAngDeg),
1825 &(opts.MinL),
1826 &(opts.NmNsFr),&(opts.Nseed),
1827 &(opts.Nmonte))) {
1828 ERROR_message("Failed to get options");
1829 exit(1);
1830 }
1831 NI_free_element(nel); nel=NULL;
1832 }
1833 }
1834
1835 if (DEF_DTI && !(opts.dti_inpref || opts.dti_listname) ){
1836 ERROR_exit("Looks like DTI, but -dti_in was NULL.\n") ;
1837 exit(5);
1838 }
1839
1840 if ( opts.dti_inpref && opts.dti_listname ) {
1841 ERROR_exit("Pick using either '-dti_listname' *or* '-dti_in'.\n") ;
1842 exit(5);
1843 }
1844
1845 if ( !DEF_DTI && ( !opts.hardi_dir || !opts.hardi_gfa ) ) {
1846 ERROR_exit("Looks like HARDI, but one of -hardi_* was NULL.\n") ;
1847 exit(5);
1848 }
1849
1850 if ( (opts.dti_inpref || opts.dti_listname) &&
1851 ( opts.hardi_dir || opts.hardi_pars || opts.hardi_gfa )) {
1852 ERROR_exit("Looks like a combo of -dti_* and -hardi_* was used.\n") ;
1853 exit(5);
1854 }
1855
1856 if ( EXTRAFILE && !DEF_DTI ) {
1857 ERROR_exit("ERROR: no extra file allowed with HARDI.\t"
1858 "Did you mean to use '-hardi_pars'?\n") ;
1859 exit(6);
1860 }
1861
1862 if( TR_MODE<0) // minip
1863 ERROR_exit("!what, no mode? Need to select from: '-mode {%s | %s | %s }'",
1864 list_modes[0],list_modes[1],list_modes[2]);
1865
1866 if( TR_MODE==1) // minip
1867 if( opts.MINI_PROB_NM<1 )
1868 ERROR_exit("in '-mode %s', one needs '-mini_num D' with D>0.\n",
1869 list_modes[TR_MODE]);
1870 else
1871 INFO_message("Running with %d mini-prob iterations",
1872 opts.MINI_PROB_NM);
1873
1874 if( (TR_MODE==1) || (TR_MODE==2) ) // minip
1875 if(!insetUC)
1876 ERROR_exit("ERROR: in '-mode %s', one needs '-uncert FILE' as well.\n",
1877 list_modes[TR_MODE]);
1878
1879 if( (TR_MODE==0) && insetUC )
1880 ERROR_exit("ERROR: in '-mode %s', one doesn't use '-uncert FILE'.\n"
1881 "\t--> Do you mean to be using '-mode %s' or '-mode %s', "
1882 "which require uncertainty?\n",
1883 list_modes[TR_MODE],list_modes[1],list_modes[2]);
1884
1885 if( (TR_MODE==2) && opts.MINI_PROB_NM )
1886 ERROR_exit("ERROR: in '-mode %s', don't assign '-mini_num' values.\n",
1887 list_modes[TR_MODE]);
1888
1889 if( (TR_MODE==2) && (opts.LOG_TYPE>=0) )
1890 ERROR_exit("ERROR: in '-mode %s', one doesn't assign '-logic' values.\n",
1891 list_modes[TR_MODE]);
1892
1893 // ** * ** * ** * ** * figure out size of vec and scal files * ** * ** *
1894
1895 // for both DTI/HARDI
1896 wild_names = (char **)calloc(N_DTI_MAX_PARS, sizeof(char *));
1897 for (j=0; j<N_DTI_MAX_PARS; ++j)
1898 wild_names[j] = (char *)calloc(32, sizeof(char));
1899 if( wild_names == NULL ){
1900 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
1901 exit(123);
1902 }
1903
1904 if(DEF_DTI) { // much reworked and cleaned up (believe it or not), Oct2016
1905
1906 insetVECS = (THD_3dim_dataset **)calloc(N_DTI_VECT,
1907 sizeof(THD_3dim_dataset *));
1908 // allocate up to max, and free extras after
1909 insetPARS = (THD_3dim_dataset **)calloc(N_DTI_MAX_PARS,
1910 sizeof(THD_3dim_dataset *));
1911
1912 if( (insetPARS == NULL) || (insetVECS == NULL) ) {
1913 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
1914 exit(123);
1915 }
1916
1917 // START of listfile set of inputs:
1918 if( opts.dti_listname ) {
1919
1920 // get both scalars (all) and vecs
1921 list_for_DTI( opts.dti_listname,
1922 insetPARS,
1923 insetVECS,
1924 &EXTRAFILE, &PARS_TOP,
1925 wild_names,
1926 1);
1927 if ( EXTRAFILE )
1928 PARS_BOT = 0;
1929
1930 } // done with list-only stuff
1931 else { // GLOB for file names, because we don't know
1932 // endings/types
1933
1934 // vecs
1935 i = glob_for_DTI_vec( opts.dti_inpref,
1936 insetVECS,
1937 hardi_pref_len);
1938 // scals
1939 i = glob_for_DTI_trac( opts.dti_inpref,
1940 insetPARS,
1941 wild_names,
1942 hardi_pref_len,
1943 &PARS_TOP,
1944 opts.NO_NONDTI_SEARCH );
1945
1946 if( opts.in_EXTRA ) {
1947 EXTRAFILE = 1;
1948 PARS_BOT = 0;
1949
1950 insetPARS[0] = THD_open_dataset(opts.in_EXTRA);
1951 if( insetPARS[0] == NULL )
1952 ERROR_exit("Can't open 'extra' listed dataset '%s': ",
1953 opts.in_EXTRA);
1954
1955 DSET_load(insetPARS[0]);
1956 CHECK_LOAD_ERROR(insetPARS[0]);
1957 fprintf(stderr,"\tFound 'extra' file '%s' to be labeled '%s'\n",
1958 opts.in_EXTRA, DTI_XTRA_LABS[0]);
1959 snprintf(wild_names[0],31,"%s", DTI_XTRA_LABS[0]);
1960 }
1961
1962 } // end of DTI glob-only stuff
1963
1964 // clean up unused
1965 for( i=PARS_TOP ; i<N_DTI_MAX_PARS ; i++ )
1966 if( insetPARS[i] ) {
1967 DSET_delete(insetPARS[i]);
1968 insetPARS[i] = NULL;
1969 free(insetPARS[i]);
1970 }
1971
1972 PARS_N = PARS_TOP - 1 + EXTRAFILE;
1973 // number of output matrices in GRID file: first two parts
1974 // are simple counts of things, and the third part is '2x'
1975 // because of having both mean and stdev for those scalars
1976 Noutmat = N_DEF_PAR_LABS + opts.EXTRA_TR_PAR +
1977 2*(PARS_N);
1978 // ParLab has length Noutmat
1979 // param_grid has 2*PARS_N+1 = Noutmat-2 blocks???????? pt,Oct,2016
1980
1981 if( EXTRAFILE ) {
1982 if( !USER_UNC_FA )
1983 unc_minfa_std*= THD_subbrick_max( insetPARS[0] , 0, 1);
1984 }
1985
1986 INFO_message("Done with scalar search, found: %d parameters\n"
1987 "\t--> so will have %d output data matrices.",
1988 PARS_N, Noutmat);
1989
1990
1991 } // end of DTI
1992 else { // HARDI section: not updated in the big Nov,2016 stuff. Perhaps later.
1993
1994 wsort = 0;
1995 nsort = 0;
1996
1997 // GLOB for file names, because we don't know endings/types, nor
1998 // number for HARDI
1999 if(opts.hardi_pars) {
2000 sprintf(tprefix,"%s*",opts.hardi_pars);
2001
2002 // this island of coding, globbing and sorting due to ZSS;
2003 // see apsearch.c program for original.
2004 wild_list = SUMA_append_replace_string(wild_list, tprefix, " ", 1);
2005 fprintf(stderr,"\n++ SEARCHING for HARDI parts with: %s",wild_list);
2006 fprintf(stdout,"\n++ SCALAR FINDINGS:\n\t");
2007 MCW_wildcards(wild_list, &nglob, &wglob );
2008 if(nglob)
2009 if( !(wsort = unique_str(wglob, nglob, wild_ci, wild_noext,
2010 &nsort, &isrt)) ) {
2011 ERROR_message("Failed to sort -hard_pars");
2012 SUMA_ifree(wild_list);
2013 MCW_free_wildcards( nglob , wglob ) ;
2014 exit(1);
2015 }
2016 }
2017 else
2018 INFO_message("Not searching for additional files,"
2019 " no '-hardi_pars' given.");
2020
2021 // we know we have GFA, so enter that
2022 PARS_N = PARS_TOP = nsort + 1; // GFA + number of files
2023 insetPARS = (THD_3dim_dataset **)calloc(PARS_N,
2024 sizeof(THD_3dim_dataset *));
2025 if( insetPARS == NULL ){
2026 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
2027 exit(123);
2028 }
2029
2030 insetPARS[0] = THD_open_dataset( opts.hardi_gfa );
2031 if( insetPARS[0] == NULL )
2032 ERROR_exit("Can't open dataset '%s'", opts.hardi_gfa);
2033 DSET_load(insetPARS[0]);
2034 CHECK_LOAD_ERROR(insetPARS[0]);
2035 snprintf(wild_names[0],31,"GFA");
2036 fprintf(stdout,"\t '%s' ",wild_names[0]);
2037
2038 // set uncert of this param to be ~2% of max, unless user
2039 // expresses value
2040 if( !USER_UNC_FA )
2041 unc_minfa_std*= THD_subbrick_max( insetPARS[0] , 0, 1);
2042
2043 if (opts.hardi_pars) {
2044
2045 jj = 1;
2046 for( ii=1 ; ii<PARS_TOP ; ii++ ) {
2047
2048 dsetn = THD_open_dataset(wglob[isrt[ii-1]]);
2049 if( DSET_NVALS(dsetn) > 1 ) {// check if it's a vec or scal
2050 DSET_delete(dsetn);
2051 dsetn=NULL;
2052 PARS_N--;
2053 }
2054 else {
2055 insetPARS[jj] = dsetn;
2056
2057 // in case of underscore between pref and var name
2058 pref_offset = 0;
2059 if( *(wsort[ii-1]+hardi_pref_len) == '_')
2060 pref_offset = 1;
2061
2062 snprintf(wild_names[jj],31,"%s",
2063 wsort[ii-1]+hardi_pref_len+pref_offset);
2064 //INFO_message("FOUND another scalar '%s' (from '%s')",
2065 // wild_names[jj],wsort[ii-1]);
2066 fprintf(stdout," '%s' ",wild_names[jj]);
2067 dsetn=NULL;
2068 if( insetPARS[jj] == NULL )
2069 ERROR_exit("Can't open dataset '%s'", wglob[isrt[ii-1]]);
2070 DSET_load(insetPARS[jj]);
2071 CHECK_LOAD_ERROR(insetPARS[jj]) ;
2072 jj++;
2073 }
2074 if( jj == N_DTI_MAX_PARS ) {
2075 INFO_message("Have reached max number of allowed"
2076 " input scalars (%d).\n"
2077 "\tWill just go with the ones gotten now.\n",
2078 N_DTI_MAX_PARS) ;
2079 PARS_N = N_DTI_MAX_PARS; // pars_top reset below
2080 break;
2081 }
2082 }
2083
2084 // free any unfilled ones, reduce stored count
2085 for( i=jj ; i<PARS_TOP ; i++ ) {
2086 // free(insetPARS[i]);
2087 insetPARS[i]=NULL;
2088 }
2089
2090 PARS_TOP = PARS_N;
2091
2092 if (isrt) free(isrt); isrt = NULL;
2093 for (i=0; i<nglob; ++i) if (wsort[i]) free(wsort[i]);
2094 free(wsort); wsort = NULL;
2095 SUMA_ifree(wild_list);
2096 MCW_free_wildcards( nglob , wglob ) ;
2097
2098 }
2099
2100 fprintf(stdout,"\n");
2101
2102 // number of output matrs in GRID
2103 Noutmat = N_DEF_PAR_LABS + opts.EXTRA_TR_PAR + 2*PARS_N;
2104 // ParLab has length Noutmat
2105 // param_grid has 2*PARS_N+1 = Noutmat-2 blocks
2106 INFO_message("Done with scalar search, found: %d parameters\n"
2107 "\t--> so will have %d output data matrices.",
2108 PARS_N, Noutmat);
2109
2110
2111 //sprintf(headerHAR.scal_n[0],"%s", wild_names[0]);
2112
2113 // need to find out max number of vectors per voxel
2114 // single file
2115 dsetn = THD_open_dataset(opts.hardi_dir);
2116 if( dsetn == NULL )
2117 ERROR_exit("Can't open dataset '%s'", opts.hardi_dir);
2118 DSET_load(dsetn) ; CHECK_LOAD_ERROR(dsetn) ;
2119
2120 N_HAR = DSET_NVALS(dsetn);
2121 if( N_HAR % 3 != 0 )
2122 ERROR_exit("Number of bricks in '-hardi_dir' set must"
2123 " be a multiple of 3,"
2124 " (which %d isn't).", N_HAR );
2125 else
2126 N_HAR/= 3;
2127
2128 insetVECS = (THD_3dim_dataset **)calloc(N_HAR,
2129 sizeof(THD_3dim_dataset *) );
2130 if( insetVECS == NULL ){
2131 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
2132 exit(123);
2133 }
2134
2135 for( i=0 ; i<N_HAR ; i++ ) {
2136 insetVECS[i] = EDIT_empty_copy(insetPARS[PARS_BOT] ) ; // single brick
2137 EDIT_add_bricklist(insetVECS[i], // make it 3D vector.
2138 2, NULL , NULL , NULL );
2139 EDIT_dset_items(insetVECS[i], // should be float anyways
2140 ADN_datum_all , MRI_float ,
2141 ADN_none ) ;
2142
2143 for( j=0 ; j<3 ; j++) {
2144 temp_vec = THD_extract_float_brick(3*i+j, dsetn) ;
2145 if( MRI_FLOAT_PTR(temp_vec) == NULL) {
2146 fprintf(stderr, "\n\n Copy failure for temp_vec.\n\n");
2147 exit(124);
2148 }
2149 EDIT_substitute_brick(insetVECS[i], j,
2150 MRI_float,
2151 MRI_FLOAT_PTR(temp_vec) );
2152 mri_clear_and_free(temp_vec);
2153 }
2154 }
2155 INFO_message("FOUND vector '%s', with %d subbricks\n"
2156 "\t--> at most %d direction(s) per voxel.",
2157 opts.hardi_dir,3*N_HAR,N_HAR);
2158
2159 DSET_delete(dsetn);
2160 free(dsetn);
2161 } // end of HARDI-only section
2162
2163 // ---------------- labels for all (DTI or HARDI) -------------------------
2164
2165 ParLab = (char **)calloc(Noutmat, sizeof(char *));
2166 for (j=0; j<Noutmat; ++j)
2167 ParLab[j] = (char *)calloc(32, sizeof(char));
2168 if( ParLab == NULL ) {
2169 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
2170 exit(121);
2171 }
2172
2173 // start loading up param labels: auto-calc'ed ones first
2174 for ( j=0; j<N_DEF_PAR_LABS; j++ )
2175 ParLab[j] = strdup(DEF_PAR_LABS[j]);
2176 if ( opts.EXTRA_TR_PAR )
2177 for ( j=0; j<N_XTR_PAR_LABS; j++ )
2178 ParLab[N_DEF_PAR_LABS + j] = strdup(XTR_PAR_LABS[j]);
2179
2180 for( i=PARS_BOT ; i<PARS_TOP ; i++ ){
2181 // 'i-PARS_BOT' goes from 0 to number-of-data-scalars-minus-1;
2182 // j: taking into account offsets and spacing for mean '%s'
2183 // and stdev 's%s'
2184 j = N_DEF_PAR_LABS + opts.EXTRA_TR_PAR + 2*(i-PARS_BOT);
2185 // mean
2186 snprintf(ParLab[j], 31, "%s", wild_names[i]);
2187 // stdev
2188 snprintf(ParLab[j+1], 31, "s%s", wild_names[i]);
2189 }
2190
2191 // ----------------- get basic info, and TV header stuff ----------------
2192
2193 if ( insetPARS[PARS_BOT] ) { // should have PARS_BOT = 0 for HARDI.
2194 Nvox = DSET_NVOX(insetPARS[PARS_BOT]) ;
2195
2196 Dim = (int *)calloc(3, sizeof(int));
2197 Dim[0] = DSET_NX(insetPARS[PARS_BOT]);
2198 Dim[1] = DSET_NY(insetPARS[PARS_BOT]);
2199 Dim[2] = DSET_NZ(insetPARS[PARS_BOT]);
2200 Ledge[0] = fabs(DSET_DX(insetPARS[PARS_BOT]));
2201 Ledge[1] = fabs(DSET_DY(insetPARS[PARS_BOT]));
2202 Ledge[2] = fabs(DSET_DZ(insetPARS[PARS_BOT]));
2203 Orig[0] = DSET_XORG(insetPARS[PARS_BOT]);
2204 Orig[1] = DSET_YORG(insetPARS[PARS_BOT]);
2205 Orig[2] = DSET_ZORG(insetPARS[PARS_BOT]);
2206 VoxVol = DSET_VOXVOL(insetPARS[PARS_BOT]);
2207
2208 // this stores the original data file orientation for later use,
2209 // as well since we convert everything to RAI temporarily, as
2210 // described below
2211 voxel_order = (char *)calloc(4, sizeof(char));
2212 TV_switch = (int *)calloc(3, sizeof(int));
2213 voxel_order[0]=ORIENT_typestr[insetPARS[PARS_BOT]->daxes->xxorient][0];
2214 voxel_order[1]=ORIENT_typestr[insetPARS[PARS_BOT]->daxes->yyorient][0];
2215 voxel_order[2]=ORIENT_typestr[insetPARS[PARS_BOT]->daxes->zzorient][0];
2216 voxel_order[3]='\0';
2217
2218 if( DEF_DTI ){
2219
2220 if ( opts.TRK_OPP_ORI ) { // [PT: June 5, 2019]
2221 headerDTI.voxel_order[0] =
2222 ORIENT_typestr[ORIENT_OPPOSITE(insetPARS[PARS_BOT]->daxes->xxorient)][0];
2223 headerDTI.voxel_order[1] =
2224 ORIENT_typestr[ORIENT_OPPOSITE(insetPARS[PARS_BOT]->daxes->yyorient)][0];
2225 headerDTI.voxel_order[2] =
2226 ORIENT_typestr[ORIENT_OPPOSITE(insetPARS[PARS_BOT]->daxes->zzorient)][0];
2227 INFO_message("Have reverse orient for any TRK files: %s", headerDTI.voxel_order );
2228 }
2229 else {
2230 headerDTI.voxel_order[0] =
2231 ORIENT_typestr[insetPARS[PARS_BOT]->daxes->xxorient][0];
2232 headerDTI.voxel_order[1] =
2233 ORIENT_typestr[insetPARS[PARS_BOT]->daxes->yyorient][0];
2234 headerDTI.voxel_order[2] =
2235 ORIENT_typestr[insetPARS[PARS_BOT]->daxes->zzorient][0];
2236 }
2237
2238 for( i=0 ; i<3 ; i++) {
2239 headerDTI.dim[i] = Dim[i];
2240 headerDTI.voxel_size[i] = Ledge[i];
2241 // will want this when outputting file later for TrackVis.
2242 TV_switch[i] = !(dset_or[i]==voxel_order[i]);
2243 }
2244 }
2245 else{ // HARDI
2246 sprintf(headerHAR.scal_n[0],"%s", wild_names[0]);
2247
2248 if ( opts.TRK_OPP_ORI ) { // [PT: June 5, 2019]
2249 headerHAR.voxel_order[0] =
2250 ORIENT_typestr[ORIENT_OPPOSITE(insetPARS[0]->daxes->xxorient)][0];
2251 headerHAR.voxel_order[1] =
2252 ORIENT_typestr[ORIENT_OPPOSITE(insetPARS[0]->daxes->yyorient)][0];
2253 headerHAR.voxel_order[2] =
2254 ORIENT_typestr[ORIENT_OPPOSITE(insetPARS[0]->daxes->zzorient)][0];
2255 INFO_message("Have reverse orient for any TRK files: %s", headerHAR.voxel_order );
2256 }
2257 else {
2258 headerHAR.voxel_order[0] =
2259 ORIENT_typestr[insetPARS[0]->daxes->xxorient][0];
2260 headerHAR.voxel_order[1] =
2261 ORIENT_typestr[insetPARS[0]->daxes->yyorient][0];
2262 headerHAR.voxel_order[2] =
2263 ORIENT_typestr[insetPARS[0]->daxes->zzorient][0];
2264 }
2265
2266 for( i=0 ; i<3 ; i++) {
2267 headerHAR.dim[i] = Dim[i];
2268 headerHAR.voxel_size[i] = Ledge[i];
2269 // will want this when outputting file later for TrackVis.
2270 TV_switch[i] = !(dset_or[i]==voxel_order[i]);
2271 }
2272
2273
2274 }
2275 }
2276
2277 // ---------------------------
2278
2279 // check tot num vox match (as proxy for dims...)
2280 if( (Nvox != nvox_rois) || ( (Nvox != nvox_unc) && !DETNET) ) {
2281 if(nvox_unc ==0)
2282 INFO_message("Are you doing something probabilistic, and\n"
2283 " did you forget to use the '-uncert' option?\n"
2284 "If not, use '-logic {AND | OR}' for"
2285 " deterministic tracking.");
2286 ERROR_exit("Input datasets do not match in number of voxels!\n"
2287 "\t'-netrois' volume has %d vox, \n"
2288 "\t'-uncert' volume has %d vox, \n"
2289 "\t'-inset' volume has %d vox, \n",nvox_rois, nvox_unc, Nvox);
2290 }
2291
2292 // Nov,2016
2293 if( ( ( voxel_order[0] != 'R' ) && ( voxel_order[0] != 'L' ) ) ||
2294 ( ( voxel_order[1] != 'A' ) && ( voxel_order[1] != 'P' ) ) ||
2295 ( ( voxel_order[2] != 'S' ) && ( voxel_order[2] != 'I' ) ) )
2296 ERROR_exit("Data must have oriention with:\n"
2297 "\tindex [0] either 'R' or 'L',\n"
2298 "\tindex [1] either 'P' or 'A', and\n"
2299 "\tindex [2] either 'I' or 'S'.\n"
2300 "\tLife is hard.");
2301
2302 // all input files must have the same orientation
2303 for( j=(PARS_BOT+1); j<PARS_TOP; ++j )
2304 i = CompareSetOrients(insetPARS[PARS_BOT], insetPARS[j]);
2305 if( DEF_DTI )
2306 for( j=0; j<3; ++j )
2307 i = CompareSetOrients(insetPARS[PARS_BOT], insetVECS[j]);
2308 else
2309 for( j=0; j<N_HAR; ++j )
2310 i = CompareSetOrients(insetPARS[PARS_BOT], insetVECS[j]);
2311
2312 if(EXTRAFILE)
2313 if( !((Dim[0] == DSET_NX(insetPARS[0])) &&
2314 (Dim[1] == DSET_NY(insetPARS[0])) &&
2315 (Dim[2] == DSET_NZ(insetPARS[0]))))
2316 ERROR_exit("Dimensions of extra set '%s' don't match those of"
2317 " the DTI prop (prefix '%s').",
2318 opts.in_EXTRA, opts.dti_inpref);
2319
2320 if( ( N_HAR>0 ) && ( N_uncert>0 ) && (N_HAR+1 != N_uncert ) )
2321 ERROR_exit("Number of `-hardi_dir's is %d, while the number of\n"
2322 " uncertainty values for the directions appears to be %d.\n"
2323 " These must match.", N_HAR, N_uncert);
2324
2325
2326 if(nvox_unc && DETNET && ( opts.MINI_PROB_NM == 0 ) ){ // @)
2327 ERROR_message("Don't -detnet AND -uncert together without -mini_num.\n");
2328 exit(1);
2329 }
2330
2331 // if(!nvox_unc && ( MINI_PROB_NM ) ){ // @)
2332 // ERROR_message("If using -mini_num, then you need -uncert as well.\n");
2333 // exit(1);
2334 //}
2335
2336 if( DETNET && (opts.LOG_TYPE==-1) ) {
2337 ERROR_message("You forgot to provide a type of logic: '-logic {AND|OR}.\n");
2338 exit(1);
2339 }
2340
2341 if(DETNET && (opts.LOG_TYPE==0) && opts.ONLY_BT){
2342 INFO_message("With '-logic OR', the '-cut_at_rois' option will"
2343 " be automically turned off (-> have '-uncut_at_rois').");
2344 opts.ONLY_BT=0;
2345 }
2346
2347 if(opts.ROIS_OUT) {
2348 for( k=0 ; k<N_nets ; k++ ) { // each netw gets own file
2349 sprintf(OUT_rois,"%s_%03d.roi.labs",opts.prefix,k);
2350 if( (fout1 = fopen(OUT_rois, "w")) == NULL) {
2351 fprintf(stderr, "Error opening file %s.",OUT_rois);
2352 exit(19);
2353 }
2354 for( i=1 ; i<=NROI[k] ; i++ ) {
2355 fprintf(fout1,"%d\t\t%d\t\t%d\n",ROI_LABELS[k][i],i,
2356 (int) pow(2,i));
2357 }
2358 fclose(fout1);
2359 }
2360 }
2361
2362 flat_matr = (float ***) calloc( N_nets, sizeof(float **) );
2363 for ( i = 0 ; i < N_nets ; i++ )
2364 flat_matr[i] = (float **) calloc( Noutmat, sizeof(float *) );
2365 for ( i = 0 ; i < N_nets ; i++ )
2366 for ( j = 0 ; j < Noutmat ; j++ )
2367 flat_matr[i][j] = (float *) calloc( NROI[i]*NROI[i], sizeof(float));
2368
2369 gdset_roi_names = (char ***)calloc(N_nets, sizeof(char **));
2370 for (i=0; i< N_nets ; i++ ) {
2371 gdset_roi_names[i] = (char **)calloc(NROI[i], sizeof(char *));
2372 for (j=0; j<NROI[i]; ++j) {
2373 gdset_roi_names[i][j] = (char *)calloc(32, sizeof(char));
2374 //snprintf(gdset_roi_names[i][j],31,"N%03d:R%d", i, ROI_LABELS[i][j]);
2375 //if( OLD_LABEL )
2376 // snprintf(gdset_roi_names[i][j],31,"N%03d:R%d", i,
2377 // ROI_LABELS[i][j]);
2378 //else{
2379 snprintf(gdset_roi_names[i][j],31,"%s",
2380 ROI_STR_LABELS[i][j+1]);
2381 //fprintf(stderr," %s ",
2382 // ROI_STR_LABELS[i][j+1]);
2383 //}
2384 }
2385 }
2386
2387 if( (flat_matr == NULL) || ( gdset_roi_names == NULL) ) {
2388 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
2389 exit(14);
2390 }
2391
2392
2393 // convert to cos of rad value for comparisons, instead of using acos()
2394 MaxAng = cos(CONV*opts.MaxAngDeg);
2395
2396 // for temp storage array, just a multiple of longest dimension!
2397 // essentially, a buffer size per tract we're making
2398 if(Dim[0] > Dim[1])
2399 ArrMax = Dim[0] * 4;
2400 else
2401 ArrMax = Dim[1] * 4;
2402 if(4*Dim[2] > ArrMax)
2403 ArrMax = Dim[2] * 4;
2404
2405 // switch to add header-- option for now, added Sept. 2012
2406 // for use with map_TrackID to map tracks to different space
2407 if(opts.RECORD_ORIG) { // @)
2408 for( i=0 ; i<3 ; i++) {
2409 if(N_HAR)
2410 headerHAR.origin[i] = Orig[i];
2411 else
2412 headerDTI.origin[i] = Orig[i];
2413 }
2414 }
2415
2416 // @) IF DOING JUST DETNET TRACK
2417 // adjust what some params mean
2418 if (DETNET){
2419
2420 if( N_nets > FOPEN_MAX)
2421 ERROR_message("You have more networks (%d) than the allowable number\n"
2422 "of open stream (%d) allowed for your computer.\n"
2423 " You should divide up the networks below this max,\n"
2424 "or spring for a better computer.\n",N_nets,FOPEN_MAX);
2425
2426 if(opts.MINI_PROB_NM)
2427 opts.Nmonte = opts.MINI_PROB_NM; // only doing 1 it
2428 else
2429 opts.Nmonte = 1;
2430
2431 opts.Nseed = opts.SeedPerV[0]*opts.SeedPerV[1]*opts.SeedPerV[2];
2432
2433 LocSeed = calloc(opts.Nseed,sizeof(LocSeed));
2434 for(i=0 ; i<opts.Nseed ; i++)
2435 LocSeed[i] = calloc(3,sizeof(float));
2436
2437 if( LocSeed == NULL ){
2438 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
2439 exit(123);
2440 }
2441
2442 // regularly spaced seeds
2443 i = 0;
2444 for( ii=0 ; ii<opts.SeedPerV[0] ; ii++ )
2445 for( jj=0 ; jj<opts.SeedPerV[1] ; jj++ )
2446 for( kk=0 ; kk<opts.SeedPerV[2] ; kk++ ) {
2447 LocSeed[i][0] = (0.5 + (float) ii)/opts.SeedPerV[0];
2448 LocSeed[i][1] = (0.5 + (float) jj)/opts.SeedPerV[1];
2449 LocSeed[i][2] = (0.5 + (float) kk)/opts.SeedPerV[2];
2450 i++;
2451 }
2452
2453 NmNsThr = 1; // thresh of 1 for stats stuff
2454
2455 N_bund = (int *)calloc(N_nets, sizeof(int));
2456 if( N_bund == NULL) {
2457 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
2458 exit(123);
2459 }
2460 if(opts.LOG_TYPE) //just for AND logic
2461 for( i=0 ; i<N_nets ; i++)
2462 N_bund[i] = FlatUHT_Len(NROI[i]);
2463 else
2464 for( i=0 ; i<N_nets ; i++)
2465 N_bund[i] = FlatUHT_Len(NROI[i]);
2466
2467 tb = (TAYLOR_BUNDLE ***) calloc( N_nets, sizeof(TAYLOR_BUNDLE **) );
2468 for ( i = 0 ; i < N_nets ; i++ ) // halftri+diag notation!!
2469 tb[i] = (TAYLOR_BUNDLE **) calloc( N_bund[i], sizeof(TAYLOR_BUNDLE *));
2470 tt = (TAYLOR_TRACT **)calloc(N_nets, sizeof(TAYLOR_TRACT *));
2471 tnet = (TAYLOR_NETWORK **)calloc(N_nets, sizeof(TAYLOR_NETWORK *));
2472 id = (int *)calloc(N_nets, sizeof(int));
2473
2474 prefix_det = calloc( N_nets,sizeof(prefix_det));
2475 for(i=0 ; i<N_nets ; i++)
2476 prefix_det[i] = calloc( THD_MAX_NAME,sizeof(char));
2477
2478 if( opts.OUTPUT_TRK)
2479 #ifdef __DragonFly__
2480 {
2481 typedef struct __FILE_public zFILE;
2482 fout0 = (FILE **)calloc(N_nets, sizeof(zFILE));
2483 }
2484 #else
2485 fout0 = (FILE **)calloc(N_nets, sizeof(FILE));
2486 #endif
2487
2488 if( (prefix_det == NULL)
2489 || (tb == NULL) || (tt == NULL) || (id == NULL)
2490 || (tnet == NULL)
2491 || ( opts.OUTPUT_TRK && (fout0 == NULL)) ) {
2492 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
2493 exit(123);
2494 }
2495 }
2496 else {
2497
2498 // Process the options a little
2499 LocSeed = calloc(opts.Nseed,sizeof(LocSeed));
2500 for(i=0 ; i<opts.Nseed ; i++)
2501 LocSeed[i] = calloc(3,sizeof(float));
2502
2503 // initial value in this case
2504 for( k=0 ; k<opts.Nseed ; k++ )
2505 for( j=0 ; j<3 ; j++ )
2506 LocSeed[k][j] = rand()*1.0/RAND_MAX;
2507
2508 // will take stats on voxels with number of tracts >=NmNsThr
2509 NmNsThr = (int) floor(opts.NmNsFr*opts.Nseed*opts.Nmonte);
2510 // lower bound is 1, and also force to be 1 if posteriori is chosen
2511 if( (NmNsThr<1) || opts.POST )
2512 NmNsThr=1;
2513 }
2514
2515 if (opts.dump_opts) {
2516 if (DETNET)
2517 nel = NI_setTractAlgOpts_M(NULL, &(opts.MinFA), &(opts.MaxAngDeg),
2518 &(opts.MinL),
2519 (opts.SeedPerV));
2520 else
2521 nel = NI_setProbTractAlgOpts_M(NULL, &(opts.MinFA), &(opts.MaxAngDeg),
2522 &(opts.MinL),
2523 &(opts.NmNsFr),&(opts.Nseed),
2524 &(opts.Nmonte) );
2525
2526 WriteTractAlgOpts(opts.prefix, nel);
2527 NI_free_element(nel); nel=NULL;
2528 }
2529
2530
2531 if(opts.POST){
2532 if(opts.DUMP_TYPE==1) {
2533 INFO_message("You asked for '-dump_rois DUMP', but also chose\n"
2534 "\t'-posteriori', so you will get binary mask DUMP,\n"
2535 "\tas well as AFNI files of numbers of tracks/voxel.");
2536 opts.DUMP_TYPE=3;
2537 }
2538 if(opts.DUMP_TYPE== -1){
2539 INFO_message("You did NOT ask for individual dump of ROIs by using\n"
2540 "\t'-dump_rois {option}' but then you DID choose\n"
2541 "\t'-posteriori', so you will get a set of AFNI files\n"
2542 "having numbers of tracks/voxel.");
2543 opts.DUMP_TYPE=2;
2544 }
2545 }
2546
2547 // have all be RAI for processing here
2548 if(TV_switch[0] || TV_switch[1] || TV_switch[2]) {
2549
2550 for( i=PARS_BOT ; i<PARS_TOP ; i++){
2551 //fprintf(stderr,"\nswitching par: %d",i);
2552 dsetn = r_new_resam_dset( insetPARS[i], NULL, 0.0, 0.0, 0.0,
2553 dset_or, RESAM_NN_TYPE, NULL, 1, 0);
2554 DSET_delete(insetPARS[i]);
2555 insetPARS[i]=dsetn;
2556 dsetn=NULL;
2557 }
2558
2559 if( DEF_DTI )
2560 ii = 3;
2561 else
2562 ii = N_HAR;
2563
2564 for( i=0 ; i<ii ; i++) {
2565 dsetn = r_new_resam_dset(insetVECS[i], NULL, 0.0, 0.0, 0.0,
2566 dset_or, RESAM_NN_TYPE, NULL, 1, 0);
2567 DSET_delete(insetVECS[i]);
2568 insetVECS[i]=dsetn;
2569 dsetn=NULL;
2570 }
2571
2572 dsetn = r_new_resam_dset( ROI_set, NULL, 0.0, 0.0, 0.0,
2573 dset_or, RESAM_NN_TYPE, NULL, 1, 0);
2574 DSET_delete(ROI_set);
2575 ROI_set=dsetn;
2576 dsetn=NULL;
2577
2578 if( N_uncert>0 ){
2579 dsetn = r_new_resam_dset( insetUC, NULL, 0.0, 0.0, 0.0,
2580 dset_or, RESAM_NN_TYPE, NULL, 1, 0);
2581 DSET_delete(insetUC);
2582 insetUC=dsetn;
2583 dsetn=NULL;
2584 }
2585
2586 if(HAVE_MASK) {
2587 dsetn = r_new_resam_dset(MASK, NULL, 0.0, 0.0, 0.0,
2588 dset_or, RESAM_NN_TYPE, NULL, 1, 0);
2589 DSET_delete(MASK);
2590 MASK=dsetn;
2591 dsetn=NULL;
2592 }
2593 }
2594
2595 // ****************************************************************
2596 // ****************************************************************
2597 // make arrays for tracking
2598 // ****************************************************************
2599 // ****************************************************************
2600
2601 // will hold indices of all voxels with actual data, i.e., takes in
2602 // (i,j,k) coor and gives index, for efficiency of storage because
2603 // lots of zeros in (Dimx, Dimy,Dimz grid).
2604 // INDEX: for afni THD_* things
2605 // INDEX2: for coor*, NETROI
2606 INDEX = (int ***) calloc( Dim[0], sizeof(int **) );
2607 for ( i = 0 ; i < Dim[0] ; i++ )
2608 INDEX[i] = (int **) calloc( Dim[1], sizeof(int *) );
2609 for ( i = 0 ; i < Dim[0] ; i++ )
2610 for ( j = 0 ; j < Dim[1] ; j++ )
2611 INDEX[i][j] = (int *) calloc( Dim[2] , sizeof(int) );
2612 INDEX2 = (int ***) calloc( Dim[0] , sizeof(int **) );
2613 for ( i = 0 ; i < Dim[0] ; i++ )
2614 INDEX2[i] = (int **) calloc( Dim[1] , sizeof(int *) );
2615 for ( i = 0 ; i < Dim[0] ; i++ )
2616 for ( j = 0 ; j < Dim[1] ; j++ )
2617 INDEX2[i][j] = (int *) calloc( Dim[2], sizeof(int) );
2618
2619 mskd = (int ***) calloc( Dim[0], sizeof(int **) );
2620 for ( i = 0 ; i < Dim[0] ; i++ )
2621 mskd[i] = (int **) calloc( Dim[1], sizeof(int *) );
2622 for ( i = 0 ; i < Dim[0] ; i++ )
2623 for ( j = 0 ; j < Dim[1] ; j++ )
2624 mskd[i][j] = (int *) calloc( Dim[2], sizeof(int) );
2625
2626 if( (INDEX == NULL) || (INDEX2 == NULL) || (mskd == NULL) ) {
2627 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
2628 exit(121);
2629 }
2630
2631 Ndata = 0;
2632 idx = 0;
2633 // determine how many voxels of actual data there are
2634 for( k=0 ; k<Dim[2] ; k++ )
2635 for( j=0 ; j<Dim[1] ; j++ )
2636 for( i=0 ; i<Dim[0] ; i++ ) {
2637 INDEX[i][j][k] = idx;
2638 if( ( (HAVE_MASK==0) &&
2639 (THD_get_voxel(insetPARS[PARS_BOT],idx,0)>EPS_V) ) ||
2640 ( HAVE_MASK && (THD_get_voxel(MASK,idx,0)>0) ) )
2641 {
2642 mskd[i][j][k]=1;
2643 Ndata+=1;
2644 INDEX2[i][j][k] = Ndata;
2645 }
2646 else
2647 INDEX2[i][j][k] = 0;
2648 idx+=1;
2649 }
2650 // now, Ndata is number of voxels of actual data, will be size of grids.
2651
2652 if( Ndata<1) {
2653 fprintf(stderr, "\n\n Too few data voxels-- wrong scale of things?\n\n");
2654 exit(1221);
2655 }
2656 //INFO_message("Ndata: %d. Nvox: %d",Ndata,Nvox);
2657
2658 DirPerVox = ( short *)calloc((Ndata+1), sizeof(short));
2659 if( DirPerVox == NULL ) {
2660 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
2661 exit(123);
2662 }
2663
2664 if(N_HAR){ // so switch from initialized defaults
2665
2666 j = Setup_Ndir_per_vox( N_HAR, Dim, mskd, INDEX, INDEX2,
2667 insetVECS, DirPerVox);
2668
2669 if(N_uncert) {
2670 UNC = calloc( (Ndata+1),sizeof(UNC));
2671 for(i=0 ; i<=Ndata ; i++) // to have all ind be >=1
2672 UNC[i] = calloc(DirPerVox[i]+1, sizeof(float));
2673 }
2674 }
2675 else{
2676 j = Setup_Ndir_per_vox( 1, Dim, mskd, INDEX, INDEX2,
2677 insetVECS, DirPerVox);
2678
2679 if(N_uncert) {
2680 // store final delta e_{12}, delta e_{12}, bias_FA and std_FA
2681 UNC = calloc( (Ndata+1),sizeof(UNC));
2682 for(i=0 ; i<=Ndata ; i++) // to have all ind be >=1
2683 UNC[i] = calloc(4, sizeof(float));
2684 }
2685 }
2686
2687 if( (N_uncert>0) && (UNC == NULL) ) {
2688 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
2689 exit(122);
2690 }
2691
2692 if( !DETNET )
2693 INFO_message("Effective Monte iterations: %d."
2694 " Fraction threshold set: %.5f\n"
2695 "\t--> Ntrack voxel threshold: %d.",
2696 opts.Nseed*opts.Nmonte,opts.NmNsFr,NmNsThr);
2697
2698 // 3 comp of V1 and FA for each data voxel
2699 // ragged for HARDI to save mem
2700 coorded = calloc( (Ndata+1),sizeof(coorded)); // to have all ind be >=1
2701 for(i=0 ; i<=Ndata ; i++)
2702 coorded[i] = calloc(3*DirPerVox[i]+1, sizeof(float));
2703 // copy for perturb
2704 copy_coorded = calloc( (Ndata+1),sizeof(copy_coorded));
2705 for(i=0 ; i<=Ndata ; i++) // to have all ind be >=1
2706 copy_coorded[i] = calloc(3*DirPerVox[i]+1, sizeof(float));
2707
2708 MAPROI = calloc( (Ndata+1),sizeof(MAPROI)); // to have all ind be >=1
2709 for(i=0 ; i<=Ndata ; i++)
2710 MAPROI[i] = calloc(N_nets,sizeof(int));
2711
2712 NETROI = (int ***) calloc( (Ndata+1), sizeof(int **) );
2713 for ( i = 0 ; i<=Ndata ; i++ )
2714 NETROI[i] = (int **) calloc( N_nets, sizeof(int *) );
2715 for ( i=0 ; i<=Ndata ; i++ )
2716 for ( j=0 ; j<N_nets ; j++ ) // jth net has NROI[j] rois
2717 NETROI[i][j] = (int *) calloc( FlatUHT_Len(NROI[j]), sizeof(int) );
2718
2719
2720 if( (coorded == NULL) || (copy_coorded == NULL)
2721 || (NETROI == NULL) || (MAPROI == NULL)) {
2722 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
2723 exit(122);
2724 }
2725
2726 Prob_grid = (int **) calloc( N_nets, sizeof(int *));
2727 for ( i = 0 ; i < N_nets ; i++ )
2728 Prob_grid[i] = (int *) calloc( FlatUHT_Len(NROI[i]), sizeof(int));
2729
2730 Prob_grid_L = (float ***) calloc( N_nets, sizeof(float **));
2731 for ( i = 0 ; i < N_nets ; i++ )
2732 Prob_grid_L[i] = (float **) calloc( FlatUHT_Len(NROI[i]), sizeof(float *));
2733 for ( i = 0 ; i < N_nets ; i++ )
2734 for ( j = 0 ; j < FlatUHT_Len(NROI[i]) ; j++ )
2735 Prob_grid_L[i][j] = (float *) calloc( 2, sizeof(float));
2736
2737 Param_grid = (float ***) calloc( N_nets, sizeof(float **));
2738 for ( i = 0 ; i < N_nets ; i++ )
2739 Param_grid[i] = (float **) calloc( FlatUHT_Len(NROI[i]), sizeof(float *));
2740 for ( i = 0 ; i < N_nets ; i++ )
2741 for ( j = 0 ; j < FlatUHT_Len(NROI[i]) ; j++ )
2742 Param_grid[i][j] = (float *) calloc( Noutmat-2, sizeof(float));
2743
2744 temp_list = ( int *)calloc((MAXNROI+1), sizeof( int));
2745 list_rois = ( int *)calloc((MAXNROI+1), sizeof( int));
2746
2747 Tforw = calloc(ArrMax,sizeof(Tforw));
2748 for(i=0 ; i<ArrMax ; i++)
2749 Tforw[i] = calloc(3,sizeof(int));
2750 Ttot = calloc(2*ArrMax,sizeof(Ttot));
2751 for(i=0 ; i<2*ArrMax ; i++)
2752 Ttot[i] = calloc(3,sizeof(int));
2753 Tback = calloc(ArrMax,sizeof(Tback));
2754 for(i=0 ; i<ArrMax ; i++)
2755 Tback[i] = calloc(3,sizeof(int));
2756 // temp storage whilst tracking, physical loc
2757 flTforw = calloc(ArrMax,sizeof(flTforw));
2758 for(i=0 ; i<ArrMax ; i++)
2759 flTforw[i] = calloc(3,sizeof(int));
2760 flTback = calloc(ArrMax,sizeof(flTback));
2761 for(i=0 ; i<ArrMax ; i++)
2762 flTback[i] = calloc(3,sizeof(int));
2763 // Nov,2016: calc for all tracking situations, to know lengths for
2764 // all
2765 flTtot = calloc(2*ArrMax,sizeof(flTtot));
2766 for(i=0 ; i<2*ArrMax ; i++)
2767 flTtot[i] = calloc(3,sizeof(float));
2768
2769 if( (flTback == NULL) || (Tforw == NULL) || (Tback == NULL)
2770 || (flTforw == NULL) || (temp_list == NULL) || (Ttot == NULL)
2771 || (list_rois == NULL) || (Param_grid == NULL)
2772 || (Prob_grid == NULL) || (Prob_grid_L == NULL)
2773 || (flTtot == NULL) ) {
2774 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
2775 exit(12);
2776 }
2777
2778 // slightly stricter now, because TrackItP includes extra endpoint
2779 // on each half tract tested
2780 ArrMax-= 1;
2781
2782 // @))
2783 if(JUMPLES) {
2784
2785 // right now, this contains space for:
2786 // + a growing param,
2787 // + and a counter of tracks per vox
2788 TROUT = (float ***) calloc( (Ndata+1), sizeof(float **) );
2789 for ( i = 0 ; i<=Ndata ; i++ )
2790 TROUT[i] = (float **) calloc( N_nets, sizeof(float *) );
2791 for ( i=0 ; i<=Ndata ; i++ )
2792 for ( j=0 ; j<N_nets ; j++ )
2793 TROUT[i][j] = (float *) calloc( 2 , sizeof(float) );
2794
2795 if( TROUT == NULL ) {
2796 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
2797 exit(14);
2798 }
2799
2800 }
2801
2802 j = Setup_Labels_Indices_Unc_M_both( Dim, mskd, INDEX,
2803 INDEX2, UNC,
2804 coorded, copy_coorded,
2805 insetPARS[PARS_BOT],
2806 DirPerVox, N_HAR,
2807 insetVECS,
2808 insetUC,
2809 unc_minei_std, unc_minfa_std,
2810 N_nets, NROI,
2811 ROI_set, MAPROI,
2812 INV_LABELS, NETROI);
2813
2814 // can free uncert dsets
2815 if(nvox_unc){
2816 DSET_delete(insetUC);
2817 free(insetUC);
2818 }
2819
2820
2821 // Aug 2014
2822 if( opts.EXTRA_TR_PAR ) {
2823
2824 INFO_message("Running three extra track parameter scalings for each"
2825 " target pair:\n"
2826 "\tNTpTarVol :div. by average target volume;\n"
2827 "\tNTpTarSA :div. by average target surface area;\n"
2828 "\tNTpTarSAFA :div. by average target surface area\n"
2829 "\t bordering suprathreshold FA (or equi-\n"
2830 "\t valent WM proxy definition).");
2831
2832 // [Nnets][Nrois+1][2 -> Vol and Surf Area and SA near highFA]
2833 ROI_SIZES = (int ***) calloc( N_nets, sizeof(int **) );
2834 for ( i=0 ; i<N_nets ; i++ )
2835 ROI_SIZES[i] = (int **) calloc( INVROI[i]+1, sizeof(int *) );
2836 for ( i=0 ; i<N_nets ; i++ )
2837 for ( j=0 ; j<INVROI[i]+1 ; j++ )
2838 ROI_SIZES[i][j] = (int *) calloc( 3 , sizeof(int) );
2839
2840 if( ROI_SIZES == NULL ) {
2841 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
2842 exit(14);
2843 }
2844
2845 i = GetSizesOfTargets( ROI_SIZES,
2846 MAPROI, Dim,
2847 mskd,
2848 INDEX2,
2849 coorded, opts.MinFA,
2850 N_nets, NROI,
2851 INV_LABELS );
2852 /*
2853 fprintf(stderr, "\n\nTEST ROI\n\n");
2854 fprintf(stderr, "VOL\n");
2855 for ( i=0 ; i<N_nets ; i++ ) {
2856 for ( j=1 ; j<INVROI[i]+1 ; j++ )
2857 fprintf(stderr, "%7d ",ROI_SIZES[i][j][0]);
2858 fprintf(stderr, "\n");
2859 }
2860 fprintf(stderr, "\n\nSA\n");
2861 for ( i=0 ; i<N_nets ; i++ ) {
2862 for ( j=1 ; j<INVROI[i]+1 ; j++ )
2863 fprintf(stderr, "%7d ",ROI_SIZES[i][j][1]);
2864 fprintf(stderr, "\n");
2865 }
2866 fprintf(stderr, "\n\nSAFA\n");
2867 for ( i=0 ; i<N_nets ; i++ ) {
2868 for ( j=1 ; j<INVROI[i]+1 ; j++ )
2869 fprintf(stderr, "%7d ",ROI_SIZES[i][j][2]);
2870 fprintf(stderr, "\n");
2871 }
2872 */
2873 }
2874
2875 // *************************************************************
2876 // *************************************************************
2877 // Beginning of main loops
2878 // *************************************************************
2879 // *************************************************************
2880
2881 Numtract = 0;
2882
2883 if(DETNET){ // @)
2884
2885 // @))) start bundle tags
2886 // counter here: halftri+diag notation!!!
2887 for( hh=0 ; hh<N_nets ; hh++)
2888 for( i=0 ; i<N_bund[hh] ; i++)
2889 tb[hh][i] = AppCreateBundle(NULL, 0, NULL); // start bundle
2890
2891 if( opts.OUTPUT_TRK )
2892 for( hh=0 ; hh<N_nets ; hh++) {
2893 sprintf(OUT_bin,"%s_%03d.trk",opts.prefix,hh); // match brick num
2894 if( (fout0[hh] = fopen(OUT_bin, "w")) == NULL) {
2895 fprintf(stderr, "Error opening file %s.",OUT_bin);
2896 exit(16);
2897 }
2898 // all outputs have same header
2899 if(N_HAR)
2900 fwrite(&headerHAR,sizeof(tv_io_header),1,fout0[hh]);
2901 else
2902 fwrite(&headerDTI,sizeof(tv_io_header),1,fout0[hh]);
2903 }
2904
2905 if (get_tract_verb()) {
2906 INFO_message("Begin tracking...");
2907 }
2908 }
2909
2910 if(opts.TARG_SURF_STOP)
2911 INFO_message("Surf-to-surf track mode ON.");
2912 if(opts.TARG_SURF_TWIXT)
2913 INFO_message("\t... and 'twixt only.");
2914
2915 if(opts.NAMEIN_thru) {
2916 thrumask = THD_open_dataset(opts.NAMEIN_thru);
2917 if( thrumask == NULL )
2918 ERROR_exit("Can't open dataset '%s': for extra set.",
2919 opts.NAMEIN_thru);
2920 DSET_load(thrumask) ; CHECK_LOAD_ERROR(thrumask) ;
2921 if( !(Dim[0] == (int) DSET_NX(thrumask) &&
2922 Dim[1] == (int) DSET_NY(thrumask) &&
2923 Dim[2] == (int) DSET_NZ(thrumask) &&
2924 N_nets == (int) DSET_NVALS(thrumask) ) )
2925 ERROR_exit("The thru_mask '%s' doesn't match up fully with the data.\n"
2926 "\tCheck the number of bricks, spatial dimensions and "
2927 "the origin.",
2928 opts.NAMEIN_thru);
2929
2930 if(TV_switch[0] || TV_switch[1] || TV_switch[2]) {
2931
2932 dsetn = r_new_resam_dset( thrumask, NULL, 0.0, 0.0, 0.0,
2933 dset_or, RESAM_NN_TYPE, NULL, 1, 0);
2934 DSET_delete(thrumask);
2935 thrumask=dsetn;
2936 dsetn=NULL;
2937 }
2938 }
2939
2940 if(TR_MODE ==2) {
2941 ni = (int) opts.Nmonte / 10.;
2942 if (ni < 2)
2943 ni = opts.Nmonte;
2944 }
2945 else if(TR_MODE==1)
2946 ni = 1;
2947 if( ni<opts.Nmonte )
2948 fprintf(stderr,"++ Tracking progress count: start ...\n");
2949 t_start = time(NULL);
2950
2951 for (gg=0 ; gg<opts.Nmonte ; gg++) {
2952
2953 if( gg>0) {// first time through is no change
2954 // relative location of each seed within voxel for this iter
2955 for( k=0 ; k<opts.Nseed ; k++ ) // @) had just minorly rearr. here
2956 for( j=0 ; j<3 ; j++ )
2957 LocSeed[k][j] = rand()*1.0/RAND_MAX;
2958
2959 if(N_HAR)
2960 j = HARDI_Perturb( Dim, mskd, INDEX, INDEX2,
2961 UNC, coorded, copy_coorded,
2962 r, DirPerVox);
2963 else
2964 j = DTI_Perturb_M( Dim, mskd, INDEX, INDEX2,
2965 UNC, coorded, copy_coorded,
2966 r,
2967 insetVECS);
2968
2969 }
2970
2971 // this is where we start the tracking for a given data set
2972 // start of Monte Carlo loop
2973 for( k=0 ; k<Dim[2] ; k++ )
2974 for( j=0 ; j<Dim[1] ; j++ )
2975 for( i=0 ; i<Dim[0] ; i++ )
2976 if( copy_coorded[INDEX2[i][j][k]][0]>=opts.MinFA )
2977 for( kk=0 ; kk<opts.Nseed ; kk++ )
2978 for( ll=0 ; ll<DirPerVox[INDEX2[i][j][k]] ; ll++ ) {
2979
2980 in[0] = i;
2981 in[1] = j;
2982 in[2] = k;
2983
2984 for( jj=0 ; jj<3 ; jj++ )
2985 physin[jj] = ((float) in[jj]+LocSeed[kk][jj])*Ledge[jj];
2986
2987 len_forw = TrackItP_NEW_M( N_HAR, DirPerVox,
2988 ll, copy_coorded,
2989 in, physin, Ledge, Dim,
2990 opts.MinFA, MaxAng,
2991 ArrMax, Tforw,
2992 flTforw, 1, phys_forw,INDEX2);
2993
2994 in[0] = i; // reset, because it's changed in TrackIt func
2995 in[1] = j;
2996 in[2] = k;
2997
2998 for( jj=0 ; jj<3 ; jj++ )
2999 physin[jj] = ((float) in[jj]+LocSeed[kk][jj])*Ledge[jj];
3000
3001 len_back = TrackItP_NEW_M( N_HAR, DirPerVox,
3002 ll, copy_coorded,
3003 in, physin, Ledge, Dim,
3004 opts.MinFA, MaxAng,
3005 ArrMax, Tback,
3006 flTback, -1, phys_back,INDEX2);
3007 // b/c of overlap of starts; includes 2ends
3008 totlen = len_forw+len_back-1;
3009 totlen_phys = phys_forw[0] + phys_back[0];
3010
3011 if( totlen_phys >= opts.MinL ) {
3012 Numtract += 1; //keeping tally of tot num of tracts
3013
3014 // glue together for simpler notation later
3015 for( n=0 ; n<len_back ; n++) { // all of this
3016 rr = len_back-n-1; // read in backward
3017 for(m=0;m<3;m++)
3018 Ttot[rr][m] = Tback[n][m];
3019
3020 }
3021 for( n=1 ; n<len_forw ; n++) { // skip first->overlap
3022 rr = n+len_back-1; // put after
3023 for(m=0;m<3;m++)
3024 Ttot[rr][m] = Tforw[n][m];
3025 }
3026
3027 vB0 = ( Ttot[totlen-1][0]<0 ) ? totlen-2 : totlen-1;
3028 //printf("%d\t%f\t\t%d\t%d\n",totlen,totlen_phys,
3029 //Ttot[0][0],Ttot[totlen-1][0]);
3030
3031 // Nov,2016: calc lengths for all!
3032 //if(DETNET){ // for now, a sep mirroring option to above
3033 // glue together for simpler notation later
3034 for( n=0 ; n<len_back ; n++) { // all of this
3035 rr = len_back-n-1; // read in backward
3036 for(m=0;m<3;m++)
3037 flTtot[rr][m] = flTback[n][m];
3038 }
3039 for( n=1 ; n<len_forw ; n++) { // skip first->overlap
3040 rr = n+len_back-1; // put after
3041 for(m=0;m<3;m++)
3042 flTtot[rr][m] = flTforw[n][m];
3043 }
3044 //}
3045
3046 // at this point now, both 0- and trL-th index
3047 // have `test index'
3048 // locations in them; these are not `tract' locs
3049 // themselves, but
3050 // voxels to check in terms of connections.
3051 // if the locations are `bad'-- either tract went
3052 // back on itself,
3053 // or tried to walk out of brain data set array,
3054 // then a -1 is
3055 // stored there, so we check against that before
3056 // looking at that
3057 // index.
3058
3059 for( hh=0 ; hh<N_nets ; hh++) {
3060
3061 // checking for NOT masks, which map MAPROI[][]=-1
3062 // if running through a not mask, split tract up
3063 // walk through once, all the way, per network
3064
3065 // initialize for this network.
3066 // first, check if initial ends are bad--
3067 // would have negative value from TrackItP only if
3068 // some badness had happened.
3069
3070 vB = ( Ttot[0][0]<0 ) ? -1 : -2;
3071 vB1 = 0;
3072 // go until hitting a NOT or penult vox
3073 while( vB+vB1 < vB0 ){
3074 // clear list of ROIs
3075 for( n=0 ; n<=MAXNROI ; n++)
3076 list_rois[n] = temp_list[n] = 0;
3077
3078 // literal start and finishes
3079 vA0 = vB+2; // for `checking' purposes
3080 // this starts at either 1 or where we ended+2
3081 vA = ( vA0==0 ) ? 1 : vA0 ; // for writing
3082 for( n=vA0 ; n<=vB0 ; n++) {
3083 rr=INDEX2[Ttot[n][0]][Ttot[n][1]][Ttot[n][2]];
3084 if( MAPROI[rr][hh]>0 )
3085 list_rois[MAPROI[rr][hh]]=1;
3086 else if( MAPROI[rr][hh]<0 )
3087 break;
3088 }
3089
3090 // keep track of where'd we gotten,
3091 // as `last good voxel'
3092 // whether we made it all the way
3093 // through track, or ran
3094 // into a NOTvox; max vB here should be vB0 still.
3095 if( n == totlen ) {
3096 vB = totlen-2;
3097 vB1 = 1; // check extra one
3098 }
3099 else {
3100 vB = n-1;
3101 vB1 = 0;
3102 }
3103
3104 // now, for this track, record
3105 // first, write shorter list of which
3106 // ones were hit
3107 m = 0;
3108 for( n=1 ; n<=NROI[hh] ; n++)
3109 if(list_rois[n]>0 ) {
3110 // values stored are 1...NROI -> 0...NROI-1
3111 // keep track of which was hit
3112 temp_list[m] = n-1;
3113 m = m+1;
3114 }
3115
3116 if ( opts.NAMEIN_thru )
3117 if( m>0 ) {
3118 doesnt_pass = 1;
3119 for( mm=vA ; mm<=vB ; mm++) {
3120 rr = INDEX[Ttot[mm][0]][Ttot[mm][1]][Ttot[mm][2]];
3121 if ( THD_get_voxel(thrumask,rr,hh) ) {
3122 doesnt_pass = 0;
3123 break;
3124 }
3125 }
3126 if( doesnt_pass )
3127 m=0; // and this will not go into next section!
3128 }
3129
3130
3131 // let's keep track of where tracts connecting
3132 // regions go.
3133 // we'll keep stats on individ ROI tracks
3134 if( m>0) {
3135
3136 //if( (opts.ONLY_BT==0) || (m==1) ) { // uncut tracts
3137 if( (opts.LOG_TYPE==0) || (m==1) ) { // OR tracts
3138 // Nov,2016: not changing this
3139 // if{} part for surf_stop
3140 // estimates for tract length.
3141 // Would have to chop tracts, etc.
3142 // Might consider doing so in
3143 // future, if need arises. So,
3144 // length on diagonals would be
3145 // full length of any tract.
3146
3147 trL = vB-vA+1; // counter of len of flTtot
3148
3149 if(trL>0) {
3150
3151 bb = 0;
3152 idx3 = MatrInd_to_FlatUHT(temp_list[bb],temp_list[bb],NROI[hh]);
3153 phys_l = 0.;
3154 for( mm=vA ; mm<=vB ; mm++) {
3155
3156 if( mm>vA) // phys len: Nov,2016
3157 phys_l+= sqrt( pow(flTtot[mm][0]-flTtot[mm-1][0],2) +
3158 pow(flTtot[mm][1]-flTtot[mm-1][1],2) +
3159 pow(flTtot[mm][2]-flTtot[mm-1][2],2));
3160
3161 rr = INDEX2[Ttot[mm][0]][Ttot[mm][1]][Ttot[mm][2]];
3162
3163 // Don't loop-> we know m==1 here! Nov,2016
3164 NETROI[rr][hh][idx3]+=1;
3165 if(NETROI[rr][hh][idx3]==NmNsThr)
3166 ss=ScoreTrackGrid_M(Param_grid,
3167 INDEX[Ttot[mm][0]][Ttot[mm][1]][Ttot[mm][2]],
3168 hh, idx3,
3169 insetPARS, PARS_BOT, PARS_TOP);
3170
3171 }
3172 // store here in case doing PROB, not DET/MINI
3173 Prob_grid[hh][idx3]+= 1;
3174 Prob_grid_L[hh][idx3][0]+= phys_l; // bundle L -> mean
3175 Prob_grid_L[hh][idx3][1]+= phys_l*phys_l; // bundle L -> stdev
3176
3177 if(DETNET && opts.LOG_TYPE==0) { // OR logic
3178 // @) if created, it will be kept; data set is just for orient/orig
3179 tt[hh] = Create_Tract_NEW(0,trL-1, flTtot+vA, id[hh],
3180 insetPARS[PARS_BOT]);
3181 ++id[hh];
3182 //tb[hh][0] = AppCreateBundle(tb[hh][0], 1, tt[hh]);
3183 tb[hh][idx3] = AppCreateBundle(tb[hh][idx3], 1, tt[hh]);
3184 tt[hh] = Free_Tracts(tt[hh], 1);
3185
3186 if( opts.OUTPUT_TRK )
3187 ss = SimpleWriteDetNetTr_M(N_HAR, fout0[hh], INDEX,
3188 insetPARS, PARS_BOT, PARS_TOP,
3189 flTtot+vA, Ttot+vA, trL,
3190 TV_switch, Dim, Ledge);
3191 }
3192 }
3193 } // end of 'if or-logic or m==1'
3194 else{
3195 // first do diagonal/individual ones,
3196 // because now we
3197 // have options for the pairwise connectors
3198
3199 // DIAGONAL
3200 phys_l = 0.;
3201 for( mm=vA ; mm<=vB ; mm++) {
3202 rr = INDEX2[Ttot[mm][0]][Ttot[mm][1]][Ttot[mm][2]];
3203
3204 if( mm>vA) // phys len: Nov,2016
3205 phys_l+= sqrt( pow(flTtot[mm][0]-flTtot[mm-1][0],2) +
3206 pow(flTtot[mm][1]-flTtot[mm-1][1],2) +
3207 pow(flTtot[mm][2]-flTtot[mm-1][2],2));
3208
3209 for( bb=0 ; bb<m ; bb++) {
3210 idx3 = MatrInd_to_FlatUHT(temp_list[bb],temp_list[bb],NROI[hh]);
3211 NETROI[rr][hh][idx3]+=1;
3212 if(NETROI[rr][hh][idx3]==NmNsThr)
3213 ss=ScoreTrackGrid_M(Param_grid,
3214 INDEX[Ttot[mm][0]][Ttot[mm][1]][Ttot[mm][2]],
3215 hh, idx3,
3216 insetPARS, PARS_BOT, PARS_TOP);
3217 }
3218 }
3219
3220 for( bb=0 ; bb<m ; bb++) { // and add length info here
3221 idx3 = MatrInd_to_FlatUHT(temp_list[bb],temp_list[bb],NROI[hh]);
3222 Prob_grid[hh][idx3]+= 1;
3223 Prob_grid_L[hh][idx3][0]+= phys_l; // bundle L -> mean
3224 Prob_grid_L[hh][idx3][1]+= phys_l*phys_l; // bundle L -> stdev
3225 }
3226
3227 // CONNECTORS: walk through mult times
3228 // just do unique connectors (we know that m>=2 here...)
3229 for( bb=0 ; bb<m ; bb++)
3230 for( cc=bb+1 ; cc<m ; cc++) {
3231
3232 if( !opts.ONLY_BT ) {
3233 start_loc = vA0;
3234 trL = vB+vB1 - vA0 + 1;
3235 }
3236 else{
3237 // 2 switches for finding ROI, and 1 for current 'FIND'
3238 onoff[0]=0; onoff[1]=0; onoff[2]=0;
3239 BreakAddCont=0;
3240
3241 start_loc = vA0; // default start
3242 trL = 0; //counter of len of flTtot
3243 // now walk through each vox, keep testing and
3244 // evaluating at each step
3245 // Nov,2016: change to not just add voxels to WM_ROI
3246 // at each step of the way, but instead to record
3247 // final start/stop
3248 for( mm=vA0 ; mm<=vB+vB1 ; mm++) {
3249 rr = INDEX2[Ttot[mm][0]][Ttot[mm][1]][Ttot[mm][2]];
3250
3251 if( MAPROI[rr][hh]==temp_list[bb]+1 ) { // hit 1
3252 onoff[0]=1;
3253 onoff[2]=1;
3254 }
3255 else if( MAPROI[rr][hh]==temp_list[cc]+1 ){ // hit 2
3256 onoff[1]=1;
3257 onoff[2]=1;
3258 }
3259 else {// a miss, could be either in b/t or outside
3260 onoff[2]=0;
3261 }
3262 switch(onoff[0]+onoff[1]) {
3263 case 2:
3264 // for default "trimmed", two opts:
3265 if(onoff[2])
3266 BreakAddCont=1; // still in last ROI
3267 else{
3268 BreakAddCont=-1; // done
3269 }
3270 // for surf-to-surf, trump last
3271 // choice and stop at
3272 // this moment, the first time the
3273 // second ROI gets hit
3274 if( opts.TARG_SURF_STOP ) {
3275 BreakAddCont=-1; // done
3276 trL++; // and include this vox
3277 // will remove later if
3278 // *TWIXT mode is on.
3279 }
3280 break;
3281 case 1:
3282 BreakAddCont=1; // in 1st or in b/t
3283 break;
3284 default:
3285 BreakAddCont=0; // just keep walking
3286 break;
3287 }
3288
3289 // are in b/t, and not at edge;
3290 // and can't include 0th track element
3291 if( (BreakAddCont==1) && mm && (mm<totlen-1) ) {
3292
3293 // catch this on the first one
3294 if( trL==0 )
3295 start_loc=mm;
3296 else if( opts.TARG_SURF_STOP )
3297 // for surf-to-surf, reset at any point
3298 // where we are still in first
3299 // region
3300 if( onoff[2] ) {
3301 start_loc=mm;
3302 trL = 0; // reset counter
3303 }
3304 trL++;
3305 }
3306 else if(BreakAddCont==-1) {// done
3307 break;
3308 }
3309 else {// unnec cond...
3310 continue;
3311 }
3312 } // end of mm
3313
3314 if( opts.TARG_SURF_TWIXT ) {
3315 // to be *between* regions only;
3316 // we only need to pull back *one* on
3317 // each side of the target connections
3318 start_loc+=1;
3319 trL-=2;
3320 }
3321 }
3322
3323 if (trL>0) {
3324
3325 phys_l = 0.;
3326 // now put tract into action
3327 idx3 = MatrInd_to_FlatUHT(temp_list[bb],temp_list[cc],NROI[hh]);
3328 for( mm=start_loc ; mm<(start_loc+trL) ; mm++) {
3329
3330 if( mm>start_loc ) // phys len: Nov,2016
3331 phys_l+= sqrt( pow(flTtot[mm][0]-flTtot[mm-1][0],2) +
3332 pow(flTtot[mm][1]-flTtot[mm-1][1],2) +
3333 pow(flTtot[mm][2]-flTtot[mm-1][2],2));
3334
3335 // get both sides of param_grid, b/c just testing one,
3336 // and param_grid is symm
3337 rr = INDEX2[Ttot[mm][0]][Ttot[mm][1]][Ttot[mm][2]];
3338
3339 NETROI[rr][hh][idx3]+=1;
3340 if(NETROI[rr][hh][idx3]==NmNsThr) {
3341 ss=ScoreTrackGrid_M( Param_grid,
3342 INDEX[Ttot[mm][0]][Ttot[mm][1]][Ttot[mm][2]],
3343 hh, idx3,
3344 insetPARS, PARS_BOT, PARS_TOP);
3345 }
3346 }
3347
3348 Prob_grid[hh][idx3]+= 1;
3349 Prob_grid_L[hh][idx3][0]+= phys_l; // bundle L -> mean
3350 Prob_grid_L[hh][idx3][1]+= phys_l*phys_l; // bundle L -> stdev
3351
3352 if(DETNET && opts.LOG_TYPE==1) {
3353 // @) if created, it will be kept
3354 tt[hh] = Create_Tract_NEW(0, trL-1, flTtot+start_loc, id[hh],
3355 insetPARS[PARS_BOT]); ++id[hh];
3356
3357 lll = MatrInd_to_FlatUHT(temp_list[bb],temp_list[cc],NROI[hh]);
3358 // AUGCHECK: think about temp_list ordering...
3359 //lll = temp_list[cc]+temp_list[bb]*NROI[hh]; // sq matr coor
3360 //lll -= (temp_list[bb]*(temp_list[bb]+1))/2; // fix for tridiag.
3361 //if(lll>N_bund[hh])
3362 //fprintf(stderr, " tb[%d][%d]=%p\n",hh, lll, tb[hh][lll]);
3363 tb[hh][lll] = AppCreateBundle(tb[hh][lll], 1, tt[hh]);
3364 tt[hh] = Free_Tracts(tt[hh], 1);
3365
3366 if( opts.OUTPUT_TRK )
3367 ss = SimpleWriteDetNetTr_M(N_HAR, fout0[hh], INDEX,
3368 insetPARS, PARS_BOT, PARS_TOP,
3369 flTtot+start_loc, Ttot+start_loc, trL,
3370 TV_switch, Dim, Ledge);
3371
3372 }
3373 } // end of trL>0
3374 } // end of cc
3375 } // end of an else
3376 }// end of 'if m>0'
3377 }
3378 }
3379 }
3380 }
3381
3382 if( TR_MODE == 2 ) {
3383 if (gg && (gg % ni == 0) )
3384 fprintf(stderr,"\t%s %3.0f%% %s -> %.2f min\n",
3385 "[", gg *10./ni,"]", (float) difftime(time(NULL),t_start)/60.);
3386 }
3387 else if( TR_MODE ==1 )
3388 fprintf(stderr,"\t%s %2d/%-2d %s -> %.2f min\n",
3389 "[", gg+1,opts.Nmonte,"]", (float) difftime(time(NULL),t_start)/60.);
3390
3391
3392 } // end of Monte Carlo loop
3393
3394 if( TR_MODE == 2 )
3395 if( ni<opts.Nmonte ){
3396 fprintf(stderr,"\t%s %3.0f%% %s -> %.2f min\n",
3397 "[", 100.,"]", (float) difftime( time(NULL) ,t_start)/60.);
3398 }
3399 else if( TR_MODE ==1 ) {
3400 fprintf(stderr,"\t%s %2d/%-2d %s -> %.2f min\n",
3401 "[", opts.Nmonte,opts.Nmonte,"]",
3402 (float) difftime(time(NULL), t_start)/60.);
3403 }
3404
3405 if(DETNET){
3406 INFO_message("Done tracking, tidying up outputs...");
3407 for( i=0 ; i<N_nets ; i++){
3408 if( opts.OUTPUT_TRK )
3409 fclose(fout0[i]); // !important to do...
3410
3411 m = 0; // count num of tracts per net
3412
3413 // Feb,2015: threshold by bundle size; default
3414 // N_tract/bundle=1, but user can increase
3415 //
3416 for (j=0; j<NROI[i]; j++)
3417 for (k=j; k<NROI[i]; k++) {
3418 jj = j*NROI[i] + k;
3419 kk = j + k*NROI[i];
3420 lll = MatrInd_to_FlatUHT(j,k,NROI[i]); // the 1-index location
3421
3422 if ( tb[i][lll]->N_tracts >= opts.THRESH_BUNDS ) { //THRESH_BUNDS
3423 snprintf( EleNameStr, 128, "%s<->%s",
3424 ROI_STR_LABELS[i][j+1], ROI_STR_LABELS[i][k+1]);
3425
3426 m+= tb[i][lll]->N_tracts;
3427
3428 if (get_tract_verb()) {
3429 INFO_message(" From tracking, net[%d],"
3430 " bund[%d] has %d tracks.",
3431 i,lll+1, tb[i][lll]->N_tracts);
3432 Show_Taylor_Bundle(tb[i][lll], NULL, 3);
3433 }
3434
3435 tnet[i] = AppAddBundleToNetwork(tnet[i], &(tb[i][lll]),
3436 jj,kk, insetPARS[PARS_BOT],
3437 EleNameStr);
3438 }
3439 else if( tb[i][lll]->N_tracts ) { // if it's nonzero, but subthr
3440 INFO_message(" Removed net[%d]'s bundle '%s<->%s' "
3441 "(only %d tracts).",
3442 i, ROI_STR_LABELS[i][j+1],
3443 ROI_STR_LABELS[i][k+1],
3444 tb[i][lll]->N_tracts);
3445
3446 rr = ByeByeBundle( j,
3447 k,
3448 i,
3449 Prob_grid,
3450 Prob_grid_L,
3451 Param_grid,
3452 Noutmat-2, // len of Param_grid 3rd dim
3453 NETROI,
3454 Ndata+1, // len of NETROI in 1st dim
3455 NROI );
3456 }
3457
3458 Free_Bundle(tb[i][lll]); // free some here... and more below
3459 }
3460
3461 INFO_message("From tracking, net[%d] has %d tracks.", i, m);
3462 }
3463
3464 for( i=0 ; i<N_nets ; i++){
3465 sprintf(prefix_det[i],"%s_%03d",opts.prefix,i);
3466
3467 if (!Write_Network(tnet[i],prefix_det[i],mode))
3468 ERROR_message("Failed to write the network.");
3469
3470 tnet[i] = Free_Network(tnet[i]);
3471 }
3472 }
3473
3474 // **************************************************************
3475 // **************************************************************
3476 // Some outputs
3477 // **************************************************************
3478 // **************************************************************
3479
3480 if(Numtract > 0 ) {
3481
3482 // apply threshold with all output stats.
3483 // threshold determined by: having more than 1 voxel in the WM-ROI
3484 // (assuming that preeetty much always there will be either ==0 or >>1)
3485 // calc mean and stdevs of different Param_grid entries
3486 for( k=0 ; k<N_nets ; k++)
3487 for( i=0 ; i<FlatUHT_Len(NROI[k]) ; i++ )
3488 if(Param_grid[k][i][0]>1.5) { // need at least 2 vox per ROI
3489 for( m=0 ; m<PARS_N ; m++) {
3490 // means
3491 Param_grid[k][i][2*m+1]/= Param_grid[k][i][0];
3492 // stdevs
3493 Param_grid[k][i][2*(m+1)]-=
3494 Param_grid[k][i][0]*pow(Param_grid[k][i][2*m+1],2);
3495 Param_grid[k][i][2*(m+1)]/= Param_grid[k][i][0]-1;
3496 Param_grid[k][i][2*(m+1)] = sqrt(Param_grid[k][i][2*(m+1)]);
3497 }
3498 Prob_grid_L[k][i][0]/= Prob_grid[k][i]; // div by num of tracts
3499 if( Prob_grid[k][i] > 1.5 ) {// need at least 2 tr per bun
3500 Prob_grid_L[k][i][1]-= Prob_grid[k][i] * pow(Prob_grid_L[k][i][0],2);
3501 Prob_grid_L[k][i][1]/= Prob_grid[k][i] - 1;
3502 if( Prob_grid_L[k][i][1] > 0 ) // 23Nov,2016: no rounding err negs->nans
3503 Prob_grid_L[k][i][1] = sqrt(Prob_grid_L[k][i][1]);
3504 else
3505 Prob_grid_L[k][i][1] = 0.;
3506 }
3507 else{
3508 INFO_message("Network [%d]:\t [%d]th WM bundle "
3509 "has only one tract!", k, i);
3510 Prob_grid_L[k][i][1] = 0.;
3511 }
3512 }
3513 else {
3514 for( m=0 ; m<PARS_N ; m++) {
3515 Param_grid[k][i][2*m+1]=0.;
3516 Param_grid[k][i][2*(m+1)]=0.;
3517 }
3518 Prob_grid[k][i] = 0;
3519 Prob_grid_L[k][i][0] = 0;
3520 Prob_grid_L[k][i][1] = 0;
3521 }
3522
3523 for( k=0 ; k<N_nets ; k++) { // each netw gets own file
3524
3525 // print out prob grid
3526 sprintf(OUT_grid,"%s_%03d.grid",opts.prefix,k); // matches brick number
3527 if( (fout1 = fopen(OUT_grid, "w")) == NULL) {
3528 fprintf(stderr, "Error opening file %s.",OUT_grid);
3529 exit(19);
3530 }
3531
3532 fprintf(fout1,"# %d # Number of network ROIs\n",NROI[k]); // N_ROIs
3533 fprintf(fout1,
3534 "# %d # Number of grid matrices\n"
3535 ,Noutmat); // Num of matrices
3536
3537 // Sept 2014: label_table stuff
3538 //if( roi_dtable ) {
3539 fprintf(fout1, "# WITH_ROI_LABELS\n");
3540 for( i=1 ; i<NROI[k] ; i++ )
3541 fprintf(fout1," %10s \t",ROI_STR_LABELS[k][i]);
3542 fprintf(fout1," %10s\n",ROI_STR_LABELS[k][i]);
3543 //}
3544
3545 for( i=1 ; i<NROI[k] ; i++ ) // labels of ROIs
3546 fprintf(fout1," %10d \t",ROI_LABELS[k][i]); // at =NROI, -> \n
3547 fprintf(fout1," %10d\n",ROI_LABELS[k][i]);
3548
3549 fprintf(fout1,"# %s\n", ParLab[0]); // NT
3550 for( i=0 ; i<NROI[k] ; i++ ) {
3551 for( j=0 ; j<NROI[k]-1 ; j++ ) // b/c we put '\n' after last one.
3552 fprintf(fout1,"%12d\t",
3553 Prob_grid[k][MatrInd_to_FlatUHT(i,j,NROI[k])]);
3554 fprintf(fout1,"%12d\n",
3555 Prob_grid[k][MatrInd_to_FlatUHT(i,j,NROI[k])]);
3556 }
3557
3558 fprintf(fout1,"# %s\n",ParLab[1]); // fNT = frac num of tracts
3559 for( i=0 ; i<NROI[k] ; i++ ) {
3560 for( j=0 ; j<NROI[k]-1 ; j++ )
3561 fprintf(fout1,"%e\t",
3562 Prob_grid[k][MatrInd_to_FlatUHT(i,j,NROI[k])]*1.0/Numtract);
3563 fprintf(fout1,"%e\n",
3564 Prob_grid[k][MatrInd_to_FlatUHT(i,j,NROI[k])]*1.0/Numtract);
3565 }
3566
3567 fprintf(fout1,"# %s\n",ParLab[2]); // PV = phys vol
3568 for( i=0 ; i<NROI[k] ; i++ ) {
3569 for( j=0 ; j<NROI[k]-1 ; j++ )
3570 fprintf(fout1,"%e\t",
3571 Param_grid[k][MatrInd_to_FlatUHT(i,j,
3572 NROI[k])][0]*VoxVol);
3573 fprintf(fout1,"%e\n",
3574 Param_grid[k][MatrInd_to_FlatUHT(i,j,NROI[k])][0]*VoxVol);
3575 }
3576
3577 fprintf(fout1,"# %s\n",ParLab[3]); // fNV = frac num of vox
3578 for( i=0 ; i<NROI[k] ; i++ ) {
3579 for( j=0 ; j<NROI[k]-1 ; j++ )
3580 fprintf(fout1,"%e\t",
3581 Param_grid[k][MatrInd_to_FlatUHT(i,j,NROI[k])][0]/Ndata);
3582 fprintf(fout1,"%e\n",
3583 Param_grid[k][MatrInd_to_FlatUHT(i,j,NROI[k])][0]/Ndata);
3584 }
3585
3586 fprintf(fout1,"# %s\n",ParLab[4]); // NV = num of vox
3587 for( i=0 ; i<NROI[k] ; i++ ) {
3588 for( j=0 ; j<NROI[k]-1 ; j++ )
3589 fprintf(fout1,"%e\t",
3590 Param_grid[k][MatrInd_to_FlatUHT(i,j,NROI[k])][0]);
3591 fprintf(fout1,"%e\n",
3592 Param_grid[k][MatrInd_to_FlatUHT(i,j,NROI[k])][0]);
3593 }
3594
3595 fprintf(fout1,"# %s\n",ParLab[5]); // BL = mean bundle lengths
3596 for( i=0 ; i<NROI[k] ; i++ ) {
3597 for( j=0 ; j<NROI[k]-1 ; j++ )
3598 fprintf(fout1,"%e\t",
3599 Prob_grid_L[k][MatrInd_to_FlatUHT(i,j,NROI[k])][0]);
3600 fprintf(fout1,"%e\n",
3601 Prob_grid_L[k][MatrInd_to_FlatUHT(i,j,NROI[k])][0]);
3602 }
3603
3604 fprintf(fout1,"# %s\n",ParLab[6]); // sBL = stdev bundle lengths
3605 for( i=0 ; i<NROI[k] ; i++ ) {
3606 for( j=0 ; j<NROI[k]-1 ; j++ )
3607 fprintf(fout1,"%e\t",
3608 Prob_grid_L[k][MatrInd_to_FlatUHT(i,j,NROI[k])][1]);
3609 fprintf(fout1,"%e\n",
3610 Prob_grid_L[k][MatrInd_to_FlatUHT(i,j,NROI[k])][1]);
3611 }
3612
3613 if (opts.EXTRA_TR_PAR) {
3614 // NTpTarVol = NT scale by ave target vol
3615 fprintf(fout1,"# %s\n",ParLab[N_DEF_PAR_LABS+0]);
3616 for( i=0 ; i<NROI[k] ; i++ ) {
3617 for( j=0 ; j<NROI[k]-1 ; j++ )
3618 fprintf(fout1,"%e\t",
3619 Prob_grid[k][MatrInd_to_FlatUHT(i,j,NROI[k])]*
3620 2.0/(ROI_SIZES[k][i+1][0] + ROI_SIZES[k][j+1][0]));
3621 fprintf(fout1,"%e\n",
3622 Prob_grid[k][MatrInd_to_FlatUHT(i,j,NROI[k])]*
3623 2.0/(ROI_SIZES[k][i+1][0] + ROI_SIZES[k][j+1][0]));
3624 }
3625
3626 // NTpTarSA = NT scale by ave target SA
3627 fprintf(fout1,"# %s\n",ParLab[N_DEF_PAR_LABS+1]);
3628 for( i=0 ; i<NROI[k] ; i++ ) {
3629 for( j=0 ; j<NROI[k]-1 ; j++ )
3630 fprintf(fout1,"%e\t",
3631 Prob_grid[k][MatrInd_to_FlatUHT(i,j,NROI[k])]*
3632 2.0/(ROI_SIZES[k][i+1][1] + ROI_SIZES[k][j+1][1]));
3633 fprintf(fout1,"%e\n",
3634 Prob_grid[k][MatrInd_to_FlatUHT(i,j,NROI[k])]*
3635 2.0/(ROI_SIZES[k][i+1][1] + ROI_SIZES[k][j+1][1]));
3636 }
3637
3638 // NTpTarSAFA = NT scale by ave target SAFA
3639 fprintf(fout1,"# %s\n",ParLab[N_DEF_PAR_LABS+2]);
3640 for( i=0 ; i<NROI[k] ; i++ ) {
3641 for( j=0 ; j<NROI[k]-1 ; j++ )
3642 fprintf(fout1,"%e\t",
3643 Prob_grid[k][MatrInd_to_FlatUHT(i,j,NROI[k])]*
3644 2.0/(ROI_SIZES[k][i+1][2] + ROI_SIZES[k][j+1][2]));
3645 fprintf(fout1,"%e\n",
3646 Prob_grid[k][MatrInd_to_FlatUHT(i,j,NROI[k])]*
3647 2.0/(ROI_SIZES[k][i+1][2] + ROI_SIZES[k][j+1][2]));
3648 }
3649
3650 } // done extra par
3651
3652 for( m=1 ; m<Noutmat-N_DEF_PAR_LABS-opts.EXTRA_TR_PAR +1; m++) {
3653 fprintf(fout1,"# %s\n",ParLab[N_DEF_PAR_LABS +
3654 opts.EXTRA_TR_PAR + m - 1]);
3655 for( i=0 ; i<NROI[k] ; i++ ) {
3656 for( j=0 ; j<NROI[k]-1 ; j++ )
3657 fprintf(fout1,"%e\t",
3658 Param_grid[k][MatrInd_to_FlatUHT(i,j,NROI[k])][m]);
3659 fprintf(fout1,"%e\n",
3660 Param_grid[k][MatrInd_to_FlatUHT(i,j,NROI[k])][m]);
3661 }
3662 }
3663
3664 fclose(fout1);
3665
3666 // @#$ recapitulate *.grid files
3667 // NB: flatmatr is flat *full* matr; *_grids are flat *UHT* matr
3668 if( DEF_DTI) {
3669 for( i=0 ; i<NROI[k] ; i++ )
3670 for( j=0 ; j<NROI[k] ; j++ ) {
3671 idx3 = MatrInd_to_FlatUHT(i,j,NROI[k]);
3672 flat_matr[k][0][i*NROI[k]+j] = Prob_grid[k][idx3];
3673 flat_matr[k][1][i*NROI[k]+j] = Prob_grid[k][idx3]*1.0/Numtract;
3674 flat_matr[k][2][i*NROI[k]+j] = Param_grid[k][idx3][0]*VoxVol;
3675 flat_matr[k][3][i*NROI[k]+j] = Param_grid[k][idx3][0]*Ndata;
3676 flat_matr[k][4][i*NROI[k]+j] = Param_grid[k][idx3][0];
3677 flat_matr[k][5][i*NROI[k]+j] = Prob_grid_L[k][idx3][0]; // Nov,2016
3678 flat_matr[k][6][i*NROI[k]+j] = Prob_grid_L[k][idx3][1];
3679
3680 if(opts.EXTRA_TR_PAR){
3681 flat_matr[k][N_DEF_PAR_LABS+0][i*NROI[k]+j] = Prob_grid[k][idx3]*
3682 2.0/(ROI_SIZES[k][i+1][0] + ROI_SIZES[k][j+1][0]);
3683 flat_matr[k][N_DEF_PAR_LABS+1][i*NROI[k]+j] = Prob_grid[k][idx3]*
3684 2.0/(ROI_SIZES[k][i+1][1] + ROI_SIZES[k][j+1][1]);
3685 flat_matr[k][N_DEF_PAR_LABS+2][i*NROI[k]+j] = Prob_grid[k][idx3]*
3686 2.0/(ROI_SIZES[k][i+1][2] + ROI_SIZES[k][j+1][2]);
3687 }
3688
3689 for( m=1 ; m<Noutmat-N_DEF_PAR_LABS-opts.EXTRA_TR_PAR+1 ;m++)
3690 flat_matr[k][N_DEF_PAR_LABS+opts.EXTRA_TR_PAR+m-1]
3691 [i*NROI[k]+j] = Param_grid[k][idx3][m];
3692 }
3693 gset = SUMA_FloatVec_to_GDSET(flat_matr[k], Noutmat,
3694 NROI[k]*NROI[k],
3695 "full", ParLab,
3696 NULL, NULL, NULL);
3697
3698 if( xyz = THD_roi_cmass(ROI_set, k, ROI_LABELS[k]+1, NROI[k], 0)){
3699 if (!(SUMA_AddGDsetNodeListElement(gset, NULL,
3700 xyz, NULL, NULL,
3701 gdset_roi_names[k],
3702 NULL, NULL,
3703 NROI[k]))) { // jan. 2014
3704 ERROR_message("Failed to add node list");
3705 exit(1);
3706 }
3707 free(xyz);
3708 }
3709 else {
3710 ERROR_message("Failed in THD_roi_cmass"); exit(1);
3711 }
3712
3713 sprintf(OUT_gdset,"%s_%03d",opts.prefix,k);
3714 GDSET_netngrlink =
3715 Network_link(SUMA_FnameGet( OUT_gdset, "f",NULL));
3716 NI_add_to_group(gset->ngr, GDSET_netngrlink);
3717 NAME_gdset = SUMA_WriteDset_ns( OUT_gdset,
3718 gset, SUMA_ASCII_NIML, 1, 0);
3719 if (!NAME_gdset && !SUMA_IS_DSET_STDXXX_FORMAT(SUMA_ASCII_NIML)) {
3720 ERROR_message("Failed to write dataset."); exit(1);
3721 } else {
3722 if (NAME_gdset) SUMA_free(NAME_gdset); NAME_gdset = NULL;
3723 }
3724 SUMA_FreeDset(gset);
3725 gset=NULL;
3726
3727 }
3728 }
3729
3730 // in order to threshold the `overall' (`0th') map;
3731 // individual ones already have
3732 // been above, just by how they were created
3733 for( i=0 ; i<=Ndata ; i++ )
3734 for( j=0 ; j<N_nets ; j++) {
3735 hh = NROI[j]*(NROI[j]+1);
3736 hh/= 2;
3737 for( k=0 ; k<hh ; k++ )
3738 if( NETROI[i][j][k] < NmNsThr)
3739 NETROI[i][j][k]=0;
3740 }
3741
3742 // output AFNI files mapping WM
3743 INFO_message("Writing output (%s, same as your input): %s ...",
3744 voxel_order,opts.prefix);
3745
3746 // some output messaging
3747 for( k=0 ; k<N_nets ; k++) {
3748 hh = 0;
3749 for( i=0 ; i<NROI[k] ; i++ )
3750 for( j=i+1 ; j<NROI[k] ; j++ ) // Don't include diags
3751 if(Prob_grid[k][MatrInd_to_FlatUHT(i,j,NROI[k])]>0){
3752 hh+=1;
3753 }
3754 if( NROI[k]>1 )
3755 INFO_message("Number of pairwise connections (bundles) in netw[%d] = %d",
3756 k,hh);
3757 }
3758
3759 if (opts.OUT_INDIPAIR) // Nov. 2014
3760 i = WriteBasicProbFiles(N_nets, Ndata, Nvox, opts.prefix,
3761 insetPARS[PARS_BOT], TV_switch, voxel_order,
3762 NROI, NETROI, mskd, INDEX2, Dim,
3763 dsetn, argc, argv,
3764 ROI_STR_LABELS, opts.DUMP_with_LABELS,
3765 roi_dtable,
3766 ROI_LABELS, opts.PAIRPOWER,
3767 opts.NIFTI_OUT);
3768 else
3769 INFO_message("Not writing out *INDI* and *PAIR* maps.");
3770
3771 if(opts.DUMP_TYPE>=0)
3772 // Prob_grid just for counting number of outputs finally
3773 i = WriteIndivProbFiles(N_nets, Ndata, Nvox, Prob_grid,
3774 opts.prefix, insetPARS[PARS_BOT],
3775 TV_switch, voxel_order,NROI,
3776 NETROI, mskd, INDEX2, Dim,
3777 dsetn, argc, argv,
3778 Param_grid, opts.DUMP_TYPE,
3779 opts.DUMP_ORIG_LABS, ROI_LABELS, opts.POST,
3780 ROI_STR_LABELS, opts.DUMP_with_LABELS,
3781 opts.NIFTI_OUT);
3782
3783 //INFO_message("Brainwide total number of tracts found = %d",Numtract);
3784 }
3785 else{
3786 INFO_message(" No Tracts Found!!!");
3787
3788 sprintf(OUT_map,"%s.pmap",opts.prefix);
3789 if( (fout1 = fopen(OUT_map, "w")) == NULL) {
3790 fprintf(stderr, "Error opening file %s.",OUT_map);
3791 exit(19);
3792 }
3793
3794 fprintf(fout1,"0!\n");
3795 fclose(fout1);
3796
3797 sprintf(OUT_grid,"%s.grid",opts.prefix);
3798 if( (fout1 = fopen(OUT_grid, "w")) == NULL) {
3799 fprintf(stderr, "Error opening file %s.",OUT_grid);
3800 exit(19);
3801 }
3802
3803 fprintf(fout1,"0!\n");
3804 fclose(fout1);
3805 }
3806
3807 // ************************************************************
3808 // ************************************************************
3809 // Freeing
3810 // ************************************************************
3811 // ************************************************************
3812
3813 if(LabTabStr)
3814 free(LabTabStr);
3815 if(roi_dtable)
3816 free(roi_dtable);
3817
3818 if( opts.EXTRA_TR_PAR ) {
3819 for ( i=0 ; i<N_nets ; i++ )
3820 for ( j=0 ; j<INVROI[i]+1 ; j++ )
3821 free(ROI_SIZES[i][j]);
3822 for ( i=0 ; i<N_nets ; i++ )
3823 free(ROI_SIZES[i]);
3824 free(ROI_SIZES);
3825 }
3826
3827 for ( i=0 ; i<N_nets ; i++ )
3828 for ( j=0 ; j<NROI[i]+1 ; j++ )
3829 free(ROI_STR_LABELS[i][j]);
3830 for ( i=0 ; i<N_nets ; i++ )
3831 free(ROI_STR_LABELS[i]);
3832 free(ROI_STR_LABELS);
3833
3834 for( i=PARS_BOT ; i<PARS_TOP ; i++){
3835 DSET_delete(insetPARS[i]);
3836 free(insetPARS[i]);
3837 }
3838 free(insetPARS);
3839
3840 if(DEF_DTI)
3841 ii = 3;
3842 else
3843 ii = N_HAR;
3844 for( i=0 ; i<ii ; i++){
3845 DSET_delete(insetVECS[i]);
3846 free(insetVECS[i]);
3847 }
3848 free(insetVECS);
3849
3850 for( i=0 ; i<Noutmat ; i++)
3851 free(ParLab[i]);
3852 free(ParLab);
3853
3854 for( i=0 ; i<N_DTI_MAX_PARS ; i++)
3855 free(wild_names[i]);
3856 free(wild_names);
3857
3858 for ( i = 0 ; i < N_nets ; i++ ) {
3859 for (j = 0; j < NROI[i]; ++j)
3860 free(gdset_roi_names[i][j]);
3861 free(gdset_roi_names[i]);
3862 }
3863 free(gdset_roi_names);
3864
3865 for ( i = 0 ; i < N_nets ; i++ )
3866 for ( j = 0 ; j < Noutmat ; j++ )
3867 free(flat_matr[i][j]);
3868 for ( i = 0 ; i < N_nets ; i++ )
3869 free(flat_matr[i]);
3870 free(flat_matr);
3871
3872 DSET_delete(ROI_set);
3873
3874 free(ROI_set);
3875 /* free(opts.prefix);
3876 free(opts.dti_inpref);
3877 free(NAMEIN_netrois);
3878 if ( NAMEIN_uncert )
3879 free(NAMEIN_uncert);
3880 if ( NAMEIN_mask )
3881 free(NAMEIN_mask);
3882 if( in_EXTRA )
3883 free( in_EXTRA );
3884 if( opts.algopt_file_name )
3885 free( opts.algopt_file_name );
3886 if ( NAMEIN_mode )
3887 free(NAMEIN_mode);
3888 if ( NAMEIN_uncEI )
3889 free( NAMEIN_uncEI );
3890 if ( NAMEIN_uncFA )
3891 free( NAMEIN_uncFA );
3892 if ( NAMEIN_outmode )
3893 free(NAMEIN_outmode);*/
3894
3895 free(DirPerVox);
3896 if(HAVE_MASK) {
3897 DSET_delete(MASK);
3898 free(MASK);
3899 }
3900
3901 for( i=0 ; i<2*ArrMax ; i++)
3902 free(Ttot[i]);
3903 free(Ttot);
3904
3905 for( i=0 ; i<ArrMax ; i++) {
3906 free(Tforw[i]);
3907 free(Tback[i]);
3908 free(flTforw[i]);
3909 free(flTback[i]);
3910 }
3911 free(Tforw);
3912 free(Tback);
3913 free(flTforw);
3914 free(flTback);
3915
3916 for( i=0 ; i<2*ArrMax ; i++) {
3917 free(flTtot[i]);
3918 }
3919 free(flTtot);
3920
3921 if(DETNET){ // @) freeing
3922 // Feb,2015: this pointer should already have been freed above!
3923 //for( i=0 ; i<N_nets ; i++)
3924 // for ( j=0 ; j<N_bund[i] ; j++ ) // halftri+diag notation!!!
3925 // free(tb[i][j]);
3926 for( i=0 ; i<N_nets ; i++) {
3927 free(tb[i]);
3928 free(tt[i]);
3929 free(tnet[i]);
3930 }
3931 free(id);
3932 free(tb);
3933 for( i=0 ; i<N_nets ; i++) {
3934 free(prefix_det[i]);
3935 }
3936 free(prefix_det);
3937 free(N_bund);
3938
3939 }
3940
3941 if(JUMPLES){
3942 for( k=0 ; k<=Ndata ; k++)
3943 for( m=0 ; m<N_nets ; m++)
3944 free(TROUT[k][m]);
3945 for( k=0 ; k<=Ndata ; k++)
3946 free(TROUT[k]);
3947 free(TROUT);
3948
3949 }
3950
3951 for( i=0 ; i<opts.Nseed ; i++)
3952 free(LocSeed[i]);
3953 free(LocSeed);
3954
3955 for( k=0 ; k<=Ndata ; k++)
3956 for( m=0 ; m<N_nets ; m++) {
3957 free(NETROI[k][m]);
3958 }
3959 for( k=0 ; k<=Ndata ; k++) {
3960 free(coorded[k]);
3961 free(copy_coorded[k]);
3962 free(NETROI[k]);
3963 free(MAPROI[k]);
3964 if(N_uncert)
3965 free(UNC[k]);
3966 }
3967 free(coorded);
3968 free(copy_coorded);
3969 free(NETROI);
3970 free(MAPROI);
3971 if(N_uncert)
3972 free(UNC);
3973
3974 for( i=0 ; i<Dim[0] ; i++)
3975 for( j=0 ; j<Dim[1] ; j++) {
3976 free(INDEX[i][j]);
3977 free(mskd[i][j]);
3978 }
3979 for( i=0 ; i<Dim[0] ; i++) {
3980 free(INDEX[i]);
3981 free(mskd[i]);
3982 }
3983 free(INDEX);
3984 free(mskd);
3985
3986 /*if(THRU_MASK>0){
3987 for( k=0 ; k<=Ndata ; k++)
3988 free(thrumask[k]);
3989 free(thrumask);
3990 }*/
3991 if(opts.NAMEIN_thru) {
3992 DSET_delete(thrumask);
3993 free(thrumask);
3994 }
3995
3996 for( i=0 ; i<Dim[0] ; i++)
3997 for( j=0 ; j<Dim[1] ; j++)
3998 free(INDEX2[i][j]);
3999 for( i=0 ; i<Dim[0] ; i++)
4000 free(INDEX2[i]);
4001 free(INDEX2);
4002
4003 for( i=0 ; i<N_nets ; i++)
4004 for( j=0 ; j<FlatUHT_Len(NROI[i]) ; j++) {
4005 free(Param_grid[i][j]);
4006 free(Prob_grid_L[i][j]);
4007 }
4008 for( i=0 ; i<N_nets ; i++) {
4009 free(Prob_grid[i]);
4010 free(Param_grid[i]);
4011 free(Prob_grid_L[i]);
4012 }
4013 free(Prob_grid);
4014 free(Param_grid);
4015 free(Prob_grid_L);
4016
4017 free(temp_list);
4018 free(list_rois);
4019
4020 for( i=0 ; i<N_nets ; i++) {
4021 free(ROI_LABELS[i]);
4022 free(INV_LABELS[i]);
4023 }
4024 free(ROI_LABELS);
4025 free(INV_LABELS);
4026
4027 free(NROI);
4028 free(INVROI);
4029 gsl_rng_free(r); // free also
4030 free(Dim);
4031 free(TV_switch);
4032 free(voxel_order);
4033
4034 return 1;
4035 }
4036
4037
4038
4039
4040
4041 // ORIGINAL GDSET_TOY example, moved from it's place among other
4042 // command line input options
4043
4044 /* if ( strcmp(argv[iarg],"-gdset_toy") == 0 ) {
4045 float x[3]={-40, 38, 4.7};
4046 float y[3]={-27, -34, 14.4};
4047 float z[3]={26, 41, 55};
4048 // alternate coordinate specification
4049 float xyz[9]={ x[0], y[0], z[0], x[1], y[1], z[1], x[2], y[2], z[2]};
4050 float **mv=NULL;
4051 int i,j,k=0;
4052 char **labs, *NameOut;
4053 SUMA_DSET *gset=NULL;
4054 TAYLOR_TRACT *tt=NULL;
4055 TAYLOR_BUNDLE *tb=NULL;
4056 TAYLOR_NETWORK *net=NULL;
4057 NI_group *netngrlink=NULL, *netngr=NULL;
4058
4059 int onexyz = 1; // set to 1 to specify coords by xyz triplets
4060 // in demo. 0 to specify x,y, and z as separate
4061 // vectors
4062
4063 INFO_message("A demonstration for creating/writing graph dset\n"
4064 "The vast majority of this content is for creating\n"
4065 "dummy data.\n"
4066 "Non dummy sections are marked with -->\n");
4067
4068 // Create 2 matrices and some matrix labels
4069 mv = (float **)calloc(2, sizeof(float*));
4070 for (i=0; i<2; ++i) mv[i] = (float *)calloc(9, sizeof(float));
4071 labs = (char **)calloc(2, sizeof(char*));
4072 for (i=0; i<2; ++i) labs[i]=calloc(64, sizeof(char));
4073
4074 for (i=0; i<2; ++i) {
4075 for (j=0; j<9; ++j) {
4076 mv[i][j] = k++;
4077 }
4078 sprintf(labs[i],"Label Of Matrix %d", i);
4079 }
4080
4081 // --> Create a graph dataset out of these two matrices
4082 gset = SUMA_FloatVec_to_GDSET(mv, 2, 9, "full", labs,
4083 NULL, NULL, NULL);
4084
4085 // --> Now add the XYZ of each graph point, two ways are possible
4086 if (onexyz) { // for when your ROI centroid coords
4087 // and in one vector of XYZ triplets
4088 if (!(SUMA_AddGDsetNodeListElement(gset, NULL,
4089 xyz, NULL, NULL,NULL, NULL, NULL, 3))) { // jan. 2014
4090 ERROR_message("Failed to add node list");
4091 exit(1);
4092 }
4093 } else {
4094 if (!(SUMA_AddGDsetNodeListElement(gset, NULL,
4095 x, y, z, NULL, NULL,
4096 NULL, 3))) { // jan. 2014
4097 ERROR_message("Failed to add node list");
4098 exit(1);
4099 }
4100 }
4101
4102 // Create some dummy bundles representing edge 1-2=5, or 2-1=7
4103 tb = NULL; net = NULL; tt = NULL;
4104 tt = (TAYLOR_TRACT *)calloc(1, sizeof(TAYLOR_TRACT));
4105 tt->id=77; tt->N_pts3=12;
4106 tt->pts = (float *)calloc(tt->N_pts3,sizeof(float));
4107 tt->pts[0]=x[1]; tt->pts[1]=y[1]; tt->pts[2]=z[1];
4108 tt->pts[3]=22; tt->pts[4]=36; tt->pts[5]=40;
4109 tt->pts[6]=22; tt->pts[7]=33; tt->pts[8]=49;
4110 tt->pts[9]=x[2]; tt->pts[10]=y[2];tt->pts[11]=z[2];
4111 tb = AppCreateBundle(tb, 1, tt);
4112 tt = Free_Tracts(tt, 1);
4113 // put another track in
4114 tt = (TAYLOR_TRACT *)calloc(1, sizeof(TAYLOR_TRACT));
4115 tt->id=78; tt->N_pts3=12;
4116 tt->pts = (float *)calloc(tt->N_pts3,sizeof(float));
4117 tt->pts[0]=x[1]; tt->pts[1]=y[1]; tt->pts[2]=z[1];
4118 tt->pts[3]=23; tt->pts[4]=35; tt->pts[5]=42;
4119 tt->pts[6]=20; tt->pts[7]=32; tt->pts[8]=51;
4120 tt->pts[9]=x[2]; tt->pts[10]=y[2];tt->pts[11]=z[2];
4121 tb = AppCreateBundle(tb, 1, tt);
4122 tt = Free_Tracts(tt, 1);
4123 // add it to network
4124 net = AppAddBundleToNetwork(net, &tb, 5, 7, NULL);
4125 // make another one for edge 0-1=1 and 1-0=3
4126 tt = (TAYLOR_TRACT *)calloc(1, sizeof(TAYLOR_TRACT));
4127 tt->id=77; tt->N_pts3=15;
4128 tt->pts = (float *)calloc(tt->N_pts3,sizeof(float));
4129 tt->pts[0]=x[0]; tt->pts[1]=y[0]; tt->pts[2]=z[0];
4130 tt->pts[3]=5; tt->pts[4]=12; tt->pts[5]=17;
4131 tt->pts[6]=16; tt->pts[7]=13; tt->pts[8]=12;
4132 tt->pts[9]=20; tt->pts[10]=16; tt->pts[11]=16;
4133 tt->pts[12]=x[1];tt->pts[13]=y[1]; tt->pts[14]=z[1];
4134 tb = AppCreateBundle(tb, 1, tt);
4135 tt = Free_Tracts(tt, 1);
4136 // add bundle to network
4137 net = AppAddBundleToNetwork(net, &tb, 1, 3, NULL);
4138
4139 // --> Now put network in graph dset
4140 netngr = Network_2_NIgr(net, 1);
4141 NI_add_to_group(gset->ngr, netngr);
4142
4143 // --> Write the graph dataset
4144 NameOut = SUMA_WriteDset_ns ("toy", gset, SUMA_ASCII_NIML, 1, 0);
4145 if (!NameOut && !SUMA_IS_DSET_STDXXX_FORMAT(SUMA_ASCII_NIML)) {
4146 ERROR_message("Failed to write dataset."); exit(1);
4147 } else {
4148 if (NameOut) SUMA_free(NameOut); NameOut = NULL;
4149 }
4150
4151 // --> Now alternately you can leave the network outside of the
4152 // graph dataset and just put a link element to it
4153 NI_remove_from_group(gset->ngr, netngr);
4154 netngrlink = Network_link("toy.network");
4155 NI_add_to_group(gset->ngr, netngrlink);
4156 NameOut = SUMA_WriteDset_ns ("toy.link",
4157 gset, SUMA_ASCII_NIML, 1, 0);
4158 if (!NameOut && !SUMA_IS_DSET_STDXXX_FORMAT(SUMA_ASCII_NIML)) {
4159 ERROR_message("Failed to write dataset."); exit(1);
4160 } else {
4161 if (NameOut) SUMA_free(NameOut); NameOut = NULL;
4162 }
4163 // And of course you need to write the tract
4164 Write_NI_Network(netngr, "toy.network", NI_TEXT_MODE);
4165 // free netngr since it is no longer tucked into dset
4166 NI_free_element(netngr); netngr = NULL;
4167
4168 // cleanup for good manners
4169 for(i=0;i<2; ++i) free(mv[i]); free(mv);
4170 for(i=0;i<2; ++i) free(labs[i]); free(labs);
4171 Free_Network(net); net = NULL;
4172
4173 INFO_message(
4174 "All done. Demo graph dset is called toy.niml.dset."
4175 "graph dset with external network spec is toy.link.niml.dset\n"
4176 "Try:\n"
4177 " suma -gdset toy.link.niml.dset\n\n"
4178 "* Open a new controller (ctrl+n), switch states ('>') to\n"
4179 "see graph dset in complementary forms.\n"
4180 "* Open 'Surface Controller' to colorize graph data, etc.\n"
4181 "Interacting with such data is possible, but not \n"
4182 "yet documented. This will come soon.\n"
4183 );
4184 exit(0);
4185 }*/
4186
4187
4188
4189
GetSizesOfTargets(int *** ROI_SIZES,int ** MAPROI,int * Dim,int *** mskd,int *** INDEX2,float ** coorded,float minFA,int N_nets,int * NROI,int ** INV_LABELS)4190 int GetSizesOfTargets( int ***ROI_SIZES,
4191 int **MAPROI, int *Dim,
4192 int ***mskd,
4193 int ***INDEX2,
4194 float **coorded, float minFA,
4195 int N_nets, int *NROI,
4196 int **INV_LABELS) {
4197
4198 int i,j,k,idx,m;
4199 int ii,jj,kk;
4200 int val,gotbnd,gotbndsa;
4201
4202 // for each vox ...
4203 for( k=0 ; k<Dim[2] ; k++ )
4204 for( j=0 ; j<Dim[1] ; j++ )
4205 for( i=0 ; i<Dim[0] ; i++ )
4206 // ... which is in the masked region ...
4207 if(mskd[i][j][k]) {
4208 idx = INDEX2[i][j][k];
4209 // ... check each network for an ROI label ...
4210 for( m=0 ; m<N_nets ; m++ ) {
4211 // ... and if found ...
4212 val = MAPROI[idx][m];
4213
4214 if( val > 0 ) {
4215 // ... then the vol increases ...
4216 ROI_SIZES[m][ val ][0]+=1;
4217
4218 gotbnd = 1;
4219 gotbndsa = 1;
4220 // ... and we check if it's a boundary ...
4221 for( ii=i-1 ; ii<=i+1 ; ii++ )
4222 for( jj=j-1 ; jj<=j+1 ; jj++ )
4223 for( kk=k-1 ; kk<=k+1 ; kk++ )
4224 // (inside FOV)
4225 if( (ii>=0) && (ii<Dim[0]) &&
4226 (jj>=0) && (jj<Dim[1]) &&
4227 (kk>=0) && (kk<Dim[2])) {
4228 // ... like this, ...
4229 if( MAPROI[INDEX2[ii][jj][kk]][m] != val ) {
4230 if (gotbnd) {// only get one
4231 ROI_SIZES[m][ val ][1]+=1;
4232 gotbnd = 0;
4233 }
4234 // ... and even whether it is
4235 // useful enough to be near
4236 // high-FA material
4237 if ( coorded[INDEX2[ii][jj][kk]][0]>minFA ) {
4238 if (gotbndsa) {// only get one
4239 ROI_SIZES[m][ val ][2]+=1;
4240 gotbndsa = 0;
4241 break; // well, break at least 1/3 loops
4242 }
4243 }
4244 }
4245 }
4246 }
4247 }
4248 }
4249
4250 RETURN(1);
4251 }
4252
4253