1 #include "mrilib.h"
2 #include "suma_objs.h"
3 #include "gicor.h"
4 #include "TrackIO.h"
5 #include "readglob.h" // need numbers of things for DTI in/out
6
7 /*
8 [PT: July 15, 2020] from include "suma_suma.h" -> "suma_objs.h"
9
10 [PT: Aug 21, 2020]
11 + remove duplicate line: #include "suma_objs.h"
12 + remove old testing print message
13 */
14
15 static int NI_tract_type = -1;
get_NI_tract_type(void)16 int get_NI_tract_type(void) {
17 if (NI_tract_type == -1) {
18 if ((NI_tract_type =
19 NI_rowtype_define( "TAYLOR_TRACT_DATUM",
20 TAYLOR_TRACT_DATUM_NIML_DEF)) < 0) {
21 ERROR_message("Failed to define NIML tract type");
22 return(-2);
23 }
24 }
25 return(NI_tract_type);
26 }
27
28 static int tract_verb = 0;
get_tract_verb(void)29 int get_tract_verb(void) { return(tract_verb); }
set_tract_verb(int v)30 void set_tract_verb(int v) { tract_verb=v; return; }
31
32 /*!
33
34 tracts_buff (TAYLOR_TRACT *) array of N_tractsbuf tracts.
35 free tracts_buff when this
36 function returns.
37
38 grid (THD_3dim_dataset *) grid defining coordinate space.
39 Only used at initialization level
40 */
AppCreateBundle(TAYLOR_BUNDLE * tbu,int N_tractsbuf,TAYLOR_TRACT * tracts_buff)41 TAYLOR_BUNDLE *AppCreateBundle(TAYLOR_BUNDLE *tbu, int N_tractsbuf,
42 TAYLOR_TRACT *tracts_buff)
43 {
44 TAYLOR_BUNDLE *tb=NULL;
45 int nn, tinb;
46 TAYLOR_TRACT *tt=NULL;
47
48 ENTRY("AppCreateBundle");
49
50 if (!tbu) {
51 tb = (TAYLOR_BUNDLE *)calloc(1,sizeof(TAYLOR_BUNDLE));
52 tb->N_allocated = 0;
53 tb->N_tracts = 0;
54 tb->N_points_private = -1;
55 tb->tract_P0_offset_private = NULL;
56 tb->bundle_ends = NULL;
57 } else {
58 tb = tbu;
59 tb->N_points_private = -1; /* reset so that this will get recomputed
60 when Bundle_N_points gets called */
61 }
62 while (N_tractsbuf > tb->N_allocated - tb->N_tracts) {
63 tb->N_allocated += 1000;
64 tb->tracts = (TAYLOR_TRACT*)realloc(tb->tracts,
65 tb->N_allocated*sizeof(TAYLOR_TRACT));
66 tb->tract_P0_offset_private = (int *)realloc(tb->tract_P0_offset_private,
67 tb->N_allocated*sizeof(int));;
68 }
69
70 if (tracts_buff && N_tractsbuf > 0) {
71 for (nn=0; nn<N_tractsbuf; ++nn) {
72 tinb = nn+tb->N_tracts;
73 tt = tb->tracts+tinb;
74 tt->id = tracts_buff[nn].id;
75 tt->N_pts3 = tracts_buff[nn].N_pts3;
76 tt->pts = (float *)calloc(tt->N_pts3, sizeof(float));
77 if (tract_verb > 1 && nn<3) {
78 fprintf(stderr,"AppCreateBundle %d , id %d, N_pts %d, pts %p\n",
79 nn, tt->id, TRACT_NPTS(tt), tracts_buff[nn].pts);
80 }
81 memcpy(tt->pts, tracts_buff[nn].pts, tt->N_pts3*sizeof(float));
82 if (tinb == 0) tb->tract_P0_offset_private[tinb] = 0;
83 else {
84 --tt; /* get previous tract */
85 tb->tract_P0_offset_private[tinb] =
86 tb->tract_P0_offset_private[tinb-1]+tt->N_pts3/3;
87 }
88 }
89 tb->N_tracts += N_tractsbuf;
90 }
91
92 RETURN(tb);
93 }
94
95
96 /*!
97 grid (THD_3dim_dataset *) Without the grid structure, coordinates
98 will not be in RAI, but in UHU, unholy units.
99
100 this works a bit differently than originally: tracts are not
101 delivered in half anymore but instead in full strings; length of
102 tract is not just given but the start and ending locations, A and B
103 (respectively and inclusively) in the array, because NOT masks might
104 trim away part of string for *some* nets but not all; finally,
105 tracts themselves will be included a bit differently-- instead of
106 having `test-ends' awkwardly separate, they will now be included in
107 the array even though they aren't full locations.
108
109 */
Create_Tract_NEW(int ptA,int ptB,float ** pts_buff,int id,THD_3dim_dataset * grid)110 TAYLOR_TRACT *Create_Tract_NEW(int ptA, int ptB, float **pts_buff,
111 int id, THD_3dim_dataset *grid)
112 {
113 TAYLOR_TRACT *tt=NULL;
114 int kk = 0, ii=0;
115 static int nwarn=0;
116 float ORIG[3], Ledge[3];
117
118 ENTRY("Create_Tract");
119
120 if (grid) {
121 if (ORIENT_typestr[grid->daxes->xxorient][0] != 'R' ||
122 ORIENT_typestr[grid->daxes->yyorient][0] != 'A' ||
123 ORIENT_typestr[grid->daxes->zzorient][0] != 'I' ) {
124 ERROR_message("Only expecting RAI grids");
125 RETURN(NULL);
126 }
127 } else {
128 if (!nwarn) {
129 WARNING_message("No grid, coordinates in UHU\n"
130 "Further messages muted\n");
131 ++nwarn;
132 }
133 }
134 tt = (TAYLOR_TRACT *)calloc(1, sizeof(TAYLOR_TRACT));
135 if (tt == NULL) {
136 ERROR_message("Failed to allocate tract");
137 RETURN(NULL);
138 }
139 tt->id = id; tt->N_pts3 = (ptB-ptA+1)*3;
140 if (!(tt->pts = (float *)calloc(tt->N_pts3, sizeof(float)))) {
141 ERROR_message("Failed to allocate pts vector");
142 Free_Tracts(tt,1); RETURN(NULL);
143 }
144
145 ORIG[0] = DSET_XORG(grid);
146 ORIG[1] = DSET_YORG(grid);
147 ORIG[2] = DSET_ZORG(grid);
148 Ledge[0] = fabs(DSET_DX(grid));
149 Ledge[1] = fabs(DSET_DY(grid));
150 Ledge[2] = fabs(DSET_DZ(grid));
151
152 kk=0;
153 if (pts_buff) { // A and B are inclusive vals
154 for (ii=ptA; ii<=ptB; ++ii) {
155 // tt->pts[kk] = pts_buff[ii][0]+DSET_XORG(grid);++kk;
156 // tt->pts[kk] = pts_buff[ii][1]+DSET_YORG(grid);++kk;
157 // tt->pts[kk] = pts_buff[ii][2]+DSET_ZORG(grid);++kk;
158 tt->pts[kk] = pts_buff[ii][0]+ORIG[0]-0.5*Ledge[0];++kk;
159 tt->pts[kk] = pts_buff[ii][1]+ORIG[1]-0.5*Ledge[1];++kk;
160 tt->pts[kk] = pts_buff[ii][2]+ORIG[2]-0.5*Ledge[2];++kk;
161 }
162 }
163
164 RETURN(tt);
165 }
166
167
168 /*! THIS VERSION JUST KEPT AROUND FOR 3dTrackID, UNTIL THAT IS FULLY
169 REMOVED.
170
171 grid (THD_3dim_dataset *) Without the grid structure, coordinates
172 will not be in RAI, but in UHU, unholy units.
173
174 ----> can prob delete now; Nov,2016
175
176 TAYLOR_TRACT *Create_Tract(int N_ptsB, float **pts_buffB,
177 int N_ptsF, float **pts_buffF,
178 int id, THD_3dim_dataset *grid)
179 {
180 TAYLOR_TRACT *tt=NULL;
181 int kk = 0, ii=0;
182 static int nwarn=0;
183
184 ENTRY("Create_Tract");
185
186 if (grid) {
187 if (ORIENT_typestr[grid->daxes->xxorient][0] != 'R' ||
188 ORIENT_typestr[grid->daxes->yyorient][0] != 'A' ||
189 ORIENT_typestr[grid->daxes->zzorient][0] != 'I' ) {
190 ERROR_message("Only expecting RAI grids");
191 RETURN(NULL);
192 }
193 } else {
194 if (!nwarn) {
195 WARNING_message("No grid, coordinates in UHU\n"
196 "Further messages muted\n");
197 ++nwarn;
198 }
199 }
200 tt = (TAYLOR_TRACT *)calloc(1, sizeof(TAYLOR_TRACT));
201 if (tt == NULL) {
202 ERROR_message("Failed to allocate tract");
203 RETURN(NULL);
204 }
205 tt->id = id; tt->N_pts3 = (N_ptsB+N_ptsF-1)*3;
206 if (!(tt->pts = (float *)calloc(tt->N_pts3, sizeof(float)))) {
207 ERROR_message("Failed to allocate pts vector");
208 Free_Tracts(tt,1); RETURN(NULL);
209 }
210 kk=0;
211 if (pts_buffB) {
212 for (ii=(N_ptsB-1); ii>0; --ii) {
213 tt->pts[kk] = pts_buffB[ii][0]+DSET_XORG(grid);++kk;
214 tt->pts[kk] = pts_buffB[ii][1]+DSET_YORG(grid);++kk;
215 tt->pts[kk] = pts_buffB[ii][2]+DSET_ZORG(grid);++kk;
216 }
217 }
218 if (pts_buffF) {
219 for (ii=0; ii<N_ptsF; ++ii) {
220 tt->pts[kk] = pts_buffF[ii][0]+DSET_XORG(grid);++kk;
221 tt->pts[kk] = pts_buffF[ii][1]+DSET_YORG(grid);++kk;
222 tt->pts[kk] = pts_buffF[ii][2]+DSET_ZORG(grid);++kk;
223 }
224 }
225 RETURN(tt);
226 }*/
227
Free_Tracts(TAYLOR_TRACT * tt,int n)228 TAYLOR_TRACT *Free_Tracts(TAYLOR_TRACT *tt, int n)
229 {
230 int i;
231
232 ENTRY("Free_Tract");
233 if (!tt) RETURN(NULL);
234 for (i=0; i<n; ++i) {
235 if (tt[i].pts) free(tt[i].pts);
236 }
237 free(tt);
238 RETURN(NULL);
239 }
240
Free_Bundle(TAYLOR_BUNDLE * tb)241 TAYLOR_BUNDLE *Free_Bundle(TAYLOR_BUNDLE *tb)
242 {
243 ENTRY("Free_Bundle");
244
245 if (!tb) RETURN(NULL);
246 tb->tracts = Free_Tracts(tb->tracts, tb->N_tracts);
247 if (tb->tract_P0_offset_private) free(tb->tract_P0_offset_private);
248 if (tb->bundle_ends) free(tb->bundle_ends);
249 free(tb);
250 RETURN(NULL);
251 }
252
Free_Network(TAYLOR_NETWORK * net)253 TAYLOR_NETWORK *Free_Network(TAYLOR_NETWORK *net)
254 {
255 TAYLOR_BUNDLE *tb=NULL;
256 int i;
257
258 ENTRY("Free_Network");
259
260 if (!net) RETURN(NULL);
261 if (net->grid) DSET_delete(net->grid);
262 if (net->FA) DSET_delete(net->FA);
263 if (net->tbv) {
264 for (i=0; i<net->N_tbv; ++i) {
265 tb = net->tbv[i];
266 if (tb) {
267 tb->tracts = Free_Tracts(tb->tracts, tb->N_tracts);
268 free(tb);
269 }
270 net->tbv[i]=NULL;
271 }
272 free(net->tbv);
273 }
274 if (net->bundle_tags) free(net->bundle_tags);
275 if (net->bundle_alt_tags) free(net->bundle_alt_tags);
276
277 free(net);
278
279 RETURN(NULL);
280 }
281
282
Show_Taylor_Tract(TAYLOR_TRACT * tt,FILE * out,int show_maxu)283 void Show_Taylor_Tract(TAYLOR_TRACT *tt, FILE *out, int show_maxu)
284 {
285 int show_max;
286 int ii=0;
287
288 ENTRY("Show_Taylor_Tract");
289 if (!out) out = stderr;
290 if (!tt) {
291 fprintf(out,"NULL tt");
292 EXRETURN;
293 }
294 fprintf(out," track id %d, Npts=%d\n", tt->id, TRACT_NPTS(tt));
295 if (show_maxu < 0) show_max = TRACT_NPTS(tt);
296 else if (show_maxu == 0) show_max = (TRACT_NPTS(tt) < 5) ? TRACT_NPTS(tt) : 5;
297 else show_max = show_maxu;
298 for (ii=0; ii<show_max; ++ii) {
299 fprintf(out, " %f %f %f\n",
300 tt->pts[3*ii], tt->pts[3*ii+1],tt->pts[3*ii+2]);
301 }
302 EXRETURN;
303 }
304
Show_Taylor_Bundle(TAYLOR_BUNDLE * tb,FILE * out,int show_maxu)305 void Show_Taylor_Bundle(TAYLOR_BUNDLE *tb, FILE *out, int show_maxu)
306 {
307 int show_max;
308 int ii=0;
309 ENTRY("Show_Taylor_Bundle");
310 if (!out) out = stderr;
311 if (!tb) {
312 fprintf(out,"NULL tb");
313 EXRETURN;
314 }
315 fprintf(out," Bundle has %d tracts, Ends %s\n",
316 tb->N_tracts, tb->bundle_ends ? tb->bundle_ends:"NULL");
317 if ((show_maxu < 0) || (tb->N_tracts < show_maxu)) show_max = tb->N_tracts;
318 else if (show_maxu == 0) show_max = (tb->N_tracts < 5) ? tb->N_tracts : 5;
319 else show_max = show_maxu;
320
321 for (ii=0; ii<show_max; ++ii) {
322 Show_Taylor_Tract(tb->tracts+ii, out, show_maxu);
323 }
324 EXRETURN;
325 }
326
Show_Taylor_Network(TAYLOR_NETWORK * net,FILE * out,int show_maxu,int show_maxub)327 void Show_Taylor_Network(TAYLOR_NETWORK *net, FILE *out,
328 int show_maxu, int show_maxub)
329 {
330 TAYLOR_BUNDLE *tb=NULL;
331 int show_max;
332 int ii=0;
333
334 ENTRY("Show_Taylor_Network");
335 if (!out) out = stderr;
336 if (!net) {
337 fprintf(out,"NULL net");
338 EXRETURN;
339 }
340 fprintf(out," Network has %d bundles\n", net->N_tbv);
341 if (show_maxu < 0) show_max = net->N_tbv;
342 else if (show_maxu == 0) show_max = (net->N_tbv < 5) ? net->N_tbv : 5;
343 else show_max = show_maxu;
344
345 for (ii=0; ii<show_max; ++ii)
346 Show_Taylor_Bundle(net->tbv[ii], out, show_maxub);
347
348 EXRETURN;
349 }
350
Tract_Length(TAYLOR_TRACT * tt)351 float Tract_Length(TAYLOR_TRACT *tt)
352 {
353 float l = -1.0, dx, dy, dz;
354 int i, N, i13, i03;
355 if (!tt) return(l);
356 N = tt->N_pts3/3;
357 l = 0.0;
358 for (i=1; i<N; ++i) {
359 i13 = 3*i; i03 = i13-3;
360 dx = tt->pts[i13 ]-tt->pts[i03 ];
361 dy = tt->pts[i13+1]-tt->pts[i03+1];
362 dz = tt->pts[i13+2]-tt->pts[i03+2];
363 l += sqrt(dx*dx+dy*dy+dz*dz);
364 }
365 return(l);
366 }
367
Bundle_N_points(TAYLOR_BUNDLE * bun,byte recalc)368 int Bundle_N_points(TAYLOR_BUNDLE *bun, byte recalc)
369 {
370 int it, nn;
371
372 if (!bun) return(-1);
373 if (!recalc && bun->N_points_private > 0) return(bun->N_points_private);
374 for (it=0, nn=0; it<bun->N_tracts; ++it) {
375 nn += bun->tracts[it].N_pts3;
376 }
377 nn /= 3;
378 bun->N_points_private = nn;
379 return(nn);
380 }
381
Network_N_points(TAYLOR_NETWORK * net,byte recalc)382 int Network_N_points(TAYLOR_NETWORK *net, byte recalc)
383 {
384 int nb, ib=0, nn=-1, it;
385
386 if (!net) return(-1);
387 if (!recalc && net->N_points_private > 0) return(net->N_points_private);
388 nn = 0;
389 for (ib=0; ib<net->N_tbv; ++ib) {
390 if (net->tbv[ib]) {
391 nb = 0;
392 for (it=0; it<net->tbv[ib]->N_tracts; ++it) {
393 nb += net->tbv[ib]->tracts[it].N_pts3;
394 }
395 net->tbv[ib]->N_points_private = nb/3;
396 nn += nb;
397 }
398 }
399 nn /= 3;
400 net->N_points_private = nn;
401 return(nn);
402 }
403
Network_N_tracts(TAYLOR_NETWORK * net,byte recalc)404 int Network_N_tracts(TAYLOR_NETWORK *net, byte recalc)
405 {
406 int nb, ib=0, nn=-1, it;
407
408 if (!net) return(-1);
409 if (!recalc && net->N_tracts_private > 0) return(net->N_tracts_private);
410 nn = 0;
411 for (ib=0; ib<net->N_tbv; ++ib) {
412 if (net->tbv[ib]) {
413 nn += net->tbv[ib]->N_tracts;
414 }
415 }
416 net->N_tracts_private = nn;
417 return(nn);
418 }
419
Network_Max_tract_length(TAYLOR_NETWORK * net,byte recalc,int * t,int * b)420 int Network_Max_tract_length(TAYLOR_NETWORK *net, byte recalc,
421 int *t, int *b)
422 {
423 int ib, it;
424
425 if (!net) return(-1);
426 if (!recalc && net->Longest_tract_length_private > 0) {
427 if (t) *t = net->Longest_tract_index_in_bundle_private;
428 if (b) *b = net->Longest_tract_bundle_index_in_network_private;
429 return(net->Longest_tract_length_private);
430 }
431 net->Longest_tract_length_private = 0;
432 for (ib=0; ib<net->N_tbv; ++ib) {
433 for (it=0; it<net->tbv[ib]->N_tracts; ++it) {
434 if (net->tbv[ib]->tracts[it].N_pts3 >
435 net->Longest_tract_length_private) {
436 net->Longest_tract_length_private = net->tbv[ib]->tracts[it].N_pts3;
437 net->Longest_tract_index_in_bundle_private = it;
438 net->Longest_tract_bundle_index_in_network_private = ib;
439 }
440 }
441 }
442 net->Longest_tract_length_private /= 3;
443
444 if (t) *t = net->Longest_tract_index_in_bundle_private;
445 if (b) *b = net->Longest_tract_bundle_index_in_network_private;
446 return(net->Longest_tract_length_private);
447 }
448
Network_N_bundles(TAYLOR_NETWORK * net)449 int Network_N_bundles(TAYLOR_NETWORK *net)
450 {
451 if (!net) return(-1);
452 return(net->N_tbv);
453 }
454
Network_PTB_to_1P(TAYLOR_NETWORK * net,int p,int t,int b)455 int Network_PTB_to_1P(TAYLOR_NETWORK *net, int p, int t, int b)
456 {
457 int PP, it, ip, ib;
458
459 ENTRY("Network_PTB_to_1P");
460
461 if (!net || p<0 || t<0 || b<0) RETURN(-1);
462
463 #if 0
464 fprintf(stderr,"P %d, T %d, B %d, \n"
465 "N_bundlesNet=%d, N_tractsBundle=%d, N_pointsTract %d\n",
466 p, t, b, b<net->N_tbv ? net->N_tbv:-1,
467 (b<net->N_tbv) ? net->tbv[b]->N_tracts:-1,
468 (b<net->N_tbv && t<net->tbv[b]->N_tracts) ?
469 net->tbv[b]->tracts[t].N_pts3/3 : -1);
470 #endif
471
472 if (b>=net->N_tbv) RETURN(-1);
473
474 if (t>=net->tbv[b]->N_tracts) RETURN(-1);
475
476 if ((3*p)>=net->tbv[b]->tracts[t].N_pts3) RETURN(-1);
477
478 PP = 0;
479 for (ib=0; ib<b; ++ib) {
480 PP += Bundle_N_points(net->tbv[ib], 0);
481 }
482 if (!net->tbv[b]->tract_P0_offset_private) {
483 for (it=0; it<t; ++it) {
484 PP += (net->tbv[b]->tracts[it].N_pts3/3);
485 }
486 } else {
487 /* use the precomputed offsets */
488 if (t > 0) PP += net->tbv[b]->tract_P0_offset_private[t];
489 }
490 PP += p;
491
492 RETURN(PP);
493 }
494
Network_TB_to_1T(TAYLOR_NETWORK * net,int t,int b)495 int Network_TB_to_1T(TAYLOR_NETWORK *net, int t, int b)
496 {
497 int it, ib, l1=0;
498
499 ENTRY("Network_TB_to_1T");
500
501 if (!net || b<0 || t<0) RETURN(-1);
502
503 #if 0
504 fprintf(stderr,"T %d, B %d, \n"
505 "N_bundlesNet=%d, N_tractsBundle=%d, N_pointsTract %d\n",
506 t, b, b<net->N_tbv ? net->N_tbv:-1,
507 (b<net->N_tbv) ? net->tbv[b]->N_tracts:-1,
508 (b<net->N_tbv && t<net->tbv[b]->N_tracts) ?
509 net->tbv[b]->tracts[t].N_pts3/3 : -1);
510 #endif
511
512 if (b>=net->N_tbv) RETURN(-1);
513
514 if (t>=net->tbv[b]->N_tracts) RETURN(-1);
515
516 l1 = 0;
517 for (ib=0; ib<b; ++ib) {
518 l1 += net->tbv[ib]->N_tracts;
519 }
520 l1 += t;
521
522 RETURN(l1);
523 }
524
Network_1T_to_TB(TAYLOR_NETWORK * net,int TT,int * t,int * b,int * PP0,int * PP1)525 int Network_1T_to_TB(TAYLOR_NETWORK *net, int TT, int *t, int *b,
526 int *PP0, int *PP1)
527 {
528 int ib;
529
530 ENTRY("Network_1T_to_TB");
531
532 if (!net || TT < 0) RETURN(-1);
533
534 ib = 0;
535 while (ib < net->N_tbv && net->tbv[ib]->N_tracts <= TT) {
536 TT -= net->tbv[ib]->N_tracts;
537 ++ib;
538 }
539 if (ib >= net->N_tbv) RETURN(-1);
540
541 /* We are in bundle ib, tract TT within ib */
542 if (b) *b = ib; if (t) *t = TT;
543
544 /* what is the 1st and last point (indexed into the network) of that tract? */
545 if (PP0) {
546 #if 1
547 *PP0 = Network_PTB_to_1P(net, 0, TT, ib);
548 #else /* No speed up, don't bother */
549 int kb;
550 for (kb = 0, *PP0 =0; kb< ib; ++kb)
551 *PP0 += Bundle_N_points(net->tbv[ib], 0);
552 for (kb=0; kb<TT; ++kb) *PP0 += net->tbv[ib]->tracts[kb].N_pts3/3;
553 #endif
554 if (PP1) *PP1 = *PP0+(net->tbv[ib]->tracts[TT].N_pts3/3 - 1);
555 }
556 RETURN(1);
557 }
558
Network_1B_to_1P(TAYLOR_NETWORK * net,int BB,int * PP1)559 int Network_1B_to_1P(TAYLOR_NETWORK *net, int BB, int *PP1)
560 {
561 int ib, PP0;
562
563 ENTRY("Network_1B_to_1P");
564
565 if (!net || BB < 0 || BB >= net->N_tbv) RETURN(-1);
566
567 ib = PP0 = 0;
568 while (ib < BB) {
569 PP0 += Bundle_N_points(net->tbv[ib], 0);
570 ++ib;
571 }
572
573 if (PP1) *PP1 = PP0 + Bundle_N_points(net->tbv[BB], 0)-1;
574
575 RETURN(PP0);
576 }
577
578 /* Go from point index into the whole network to
579 (bundle in network, tract in bundle, point in tract)
580 and if l1u is not NULL, return the index into the whole
581 network of the tract containing the point index.
582 The latter index is useful for retrieving data stored
583 per tract, rather than per point (SUMA_LEV1_DAT). */
Network_1P_to_PTB(TAYLOR_NETWORK * net,int PP,int * p,int * t,int * b,int * l1u)584 int Network_1P_to_PTB(TAYLOR_NETWORK *net, int PP,
585 int *p, int *t, int *b, int *l1u)
586 {
587 int ib, bnp, it, tnp, l1;
588
589 ENTRY("Network_1P_to_PTB");
590
591 if (!net || PP<0) RETURN(-1);
592
593 ib = l1 = 0;
594 while (ib < net->N_tbv && ((bnp = Bundle_N_points(net->tbv[ib], 0)) <= PP)) {
595 PP -= bnp;
596 l1 += net->tbv[ib]->N_tracts;
597 ++ib;
598 }
599 if (ib >= net->N_tbv) RETURN(-1);
600
601 /* We're in bundle ib, get tract in question */
602 it = 0;
603 while (it < net->tbv[ib]->N_tracts &&
604 ((tnp = net->tbv[ib]->tracts[it].N_pts3/3) <= PP)) {
605 PP -= tnp;
606 ++it;
607 }
608 if (it >= net->tbv[ib]->N_tracts) RETURN(-1);
609 l1 += it; /* Tract index (into whole network) */
610
611 *p = PP;
612 *t = it;
613 *b = ib;
614
615 if (l1u) *l1u= l1;
616
617 RETURN(1);
618 }
619
Tract_2_NIel(TAYLOR_TRACT * tt)620 NI_element *Tract_2_NIel(TAYLOR_TRACT *tt)
621 {
622 NI_element *nel=NULL;
623 char colabs[1024]={""};
624
625 ENTRY("Tract_2_NIel");
626
627 if (!tt || TRACT_NPTS(tt) < 0) RETURN(nel);
628
629 nel = NI_new_data_element("tract", TRACT_NPTS(tt));
630 NI_SETA_INT(nel, "id", tt->id);
631
632 if (tt->pts) {
633 strncat(colabs, "x;", (1023-strlen(colabs))*sizeof(char));
634 NI_add_column_stride(nel, NI_FLOAT, tt->pts , 3);
635 strncat(colabs, "y;", (1023-strlen(colabs))*sizeof(char));
636 NI_add_column_stride(nel, NI_FLOAT, tt->pts+1, 3);
637 strncat(colabs, "z;", (1023-strlen(colabs))*sizeof(char));
638 NI_add_column_stride(nel, NI_FLOAT, tt->pts+2, 3);
639 }
640
641 NI_set_attribute(nel,"Column_Labels", colabs);
642 RETURN(nel);
643 }
644
Tracts_2_NIel(TAYLOR_TRACT * tt,int N_tt)645 NI_element *Tracts_2_NIel(TAYLOR_TRACT *tt, int N_tt)
646 {
647 NI_element *nel=NULL;
648
649 ENTRY("Tracts_2_NIel");
650
651 if (!tt || !N_tt) RETURN(nel);
652
653 nel = NI_new_data_element("tracts", N_tt);
654 NI_add_column( nel , get_NI_tract_type(), tt );
655
656 NI_set_attribute(nel,"Column_Labels", "TaylorTract");
657 RETURN(nel);
658 }
659
NIel_2_Tracts(NI_element * nel,int * N_tracts)660 TAYLOR_TRACT *NIel_2_Tracts(NI_element *nel, int *N_tracts)
661 {
662 TAYLOR_TRACT *tt = NULL, *ttn=NULL;
663 float *fv0=NULL, *fv1=NULL, *fv2=NULL;
664 int ii=0, kk=0, nn=0;
665
666 ENTRY("NIel_2_Tracts");
667
668 *N_tracts = 0;
669 if (!nel) RETURN(tt);
670
671 if (!strcmp(nel->name,"tract")) {
672 *N_tracts = 1;
673 tt = (TAYLOR_TRACT *)calloc(*N_tracts,sizeof(TAYLOR_TRACT));
674 NI_GETA_INT(nel, "id", tt->id);
675 tt->N_pts3 = 3*nel->vec_len;
676 if (nel->vec_num >= 3) {
677 if (!(tt->pts = (float*)calloc(3*nel->vec_len, sizeof(float)))) {
678 ERROR_message("Failed to allocate");
679 Free_Tracts(tt,*N_tracts); RETURN(NULL);
680 }
681 fv0 = (float*)nel->vec[0];
682 fv1 = (float*)nel->vec[1];
683 fv2 = (float*)nel->vec[2];
684 kk=0;
685 for (ii=0; ii<TRACT_NPTS(tt); ++ii) {
686 tt->pts[kk] = fv0[ii]; ++kk;
687 tt->pts[kk] = fv1[ii]; ++kk;
688 tt->pts[kk] = fv2[ii]; ++kk;
689 }
690 }
691 } else if (!strcmp(nel->name,"tracts")) {
692 if (nel->vec_typ[0] != get_NI_tract_type()) {
693 ERROR_message("Bad vec_type, have %d, expected %d",
694 nel->vec_typ[0], get_NI_tract_type());
695 RETURN(NULL);
696 }
697 *N_tracts = nel->vec_len;
698 tt = (TAYLOR_TRACT *)calloc(*N_tracts,sizeof(TAYLOR_TRACT));
699 for (nn=0; nn<*N_tracts; ++nn) {
700 ttn = (TAYLOR_TRACT *)(nel->vec[0])+nn;
701 tt[nn].id = ttn->id;
702 tt[nn].N_pts3 = ttn->N_pts3;
703 tt[nn].pts = (float *)calloc(ttn->N_pts3, sizeof(float));
704 if (tract_verb && nn<3) {
705 fprintf(stderr,"NIel_2_Tracts %d , id %d, N_pts %d, pts %p\n",
706 nn, ttn->id, TRACT_NPTS(ttn), ttn->pts);
707 }
708 memcpy(tt[nn].pts, ttn->pts, ttn->N_pts3*sizeof(float));
709 }
710 }
711 RETURN(tt);
712 }
713
Network_2_NIgr(TAYLOR_NETWORK * net,int mode)714 NI_group *Network_2_NIgr(TAYLOR_NETWORK *net, int mode)
715 {
716 NI_element *nel=NULL;
717 NI_group *ngr=NULL, *ngrgrid=NULL, *ngrfa=NULL;
718 TAYLOR_BUNDLE *tb=NULL;
719 int ii=0, N_All_tracts, ei, ei_alt, bb;
720
721 ENTRY("Network_2_NIgr");
722
723 if ( !net || !net->tbv || net->N_tbv < 1) RETURN(ngr);
724
725 ngr = NI_new_group_element(); NI_rename_group(ngr,"network");
726 for (N_All_tracts=0, bb=0; bb<net->N_tbv; ++bb) {
727 if ((tb = net->tbv[bb])) {
728 N_All_tracts += tb->N_tracts;
729 }
730 }
731 NI_SETA_INT(ngr, "N_tracts", N_All_tracts);
732 for (bb=0; bb<net->N_tbv; ++bb) {
733 if ((tb = net->tbv[bb])) {
734 if (net->bundle_tags) ei = net->bundle_tags[bb];
735 else ei = bb;
736 if (net->bundle_alt_tags) ei_alt = net->bundle_alt_tags[bb];
737 else ei_alt = -1;
738 if (tb->tracts) {
739 if (mode == 0) { /* slow, handy */
740 for (ii=0; ii<tb->N_tracts; ++ii) {
741 nel = Tract_2_NIel(tb->tracts+ii);
742 NI_add_to_group(ngr, nel);
743 }
744 } else if (mode == 1) { /* fast */
745 nel = Tracts_2_NIel(tb->tracts, tb->N_tracts);
746 NI_SET_INT(nel,"Bundle_Tag", ei);
747 if (ei_alt >= 0) NI_SET_INT(nel,"Bundle_Alt_Tag", ei_alt);
748 if (tb->bundle_ends)
749 NI_SET_STR(nel,"Bundle_Ends", tb->bundle_ends);
750 NI_add_to_group(ngr, nel);
751 }
752 }
753
754 }
755 }
756 if (net->grid) {
757 ngrgrid = THD_dataset_to_niml(net->grid);
758 NI_set_attribute(ngrgrid,"bundle_aux_dset","grid");
759 NI_add_to_group(ngr, ngrgrid);
760 if (net->atlas_space)
761 NI_set_attribute(ngr,"atlas_space", net->atlas_space);
762 }
763 if (net->FA) {
764 ngrfa = THD_dataset_to_niml(net->FA);
765 NI_set_attribute(ngrfa,"bundle_aux_dset","FA");
766 NI_add_to_group(ngr, ngrfa);
767 }
768
769
770
771 RETURN(ngr);
772 }
773
Network_link(char * filename)774 NI_group *Network_link(char *filename)
775 {
776 NI_group *ngr=NULL;
777 char *fext = NULL;
778
779 ENTRY("Network_link");
780
781 if ( !filename) RETURN(ngr);
782 fext = SUMA_Extension(filename, ".niml.tract", NOPE);
783 ngr = NI_new_group_element(); NI_rename_group(ngr,"network_link");
784 NI_set_attribute(ngr, "network_file", fext);
785 SUMA_free(fext);
786
787 RETURN(ngr);
788 }
789
790
NIgr_2_Network(NI_group * ngr)791 TAYLOR_NETWORK *NIgr_2_Network(NI_group *ngr)
792 {
793 TAYLOR_NETWORK *net=NULL;
794 TAYLOR_BUNDLE *tbb=NULL;
795 TAYLOR_TRACT *tt=NULL;
796 NI_element *nel=NULL;
797 int ip=0, N_tracts=0, ei=0;
798 char *bad=NULL, *sbuf=NULL;
799 char tb_ends[128];
800
801 ENTRY("NIgr_2_Network");
802
803 if (!ngr) RETURN(net);
804 if (!strcmp(ngr->name,"bundle") ||/* old style */
805 !strcmp(ngr->name,"network") ) {
806 net = (TAYLOR_NETWORK *)calloc(1,sizeof(TAYLOR_NETWORK));
807 net->N_points_private = -1;
808 net->N_tracts_private = -1;
809 tbb = (TAYLOR_BUNDLE *)calloc(1,sizeof(TAYLOR_BUNDLE));
810 tbb->N_points_private = -1;
811 for( ip=0 ; ip < ngr->part_num ; ip++ ){
812 switch( ngr->part_typ[ip] ){
813 case NI_GROUP_TYPE:
814 if (!(bad = NI_get_attribute(ngr,"bundle_aux_dset"))) {
815 WARNING_message("Got unknown group in here! Plodding along");
816 }
817 if (!strcmp(bad,"grid")) {
818 net->grid = THD_niml_to_dataset((NI_group*)ngr->part[ip], 0);
819 } else if (!strcmp(bad,"FA")) {
820 net->FA = THD_niml_to_dataset((NI_group*)ngr->part[ip], 0);
821 } else {
822 WARNING_message("Not ready to feel the love for %s\n", bad);
823 }
824 if ((sbuf = NI_get_attribute((NI_group*)ngr->part[ip]
825 ,"atlas_space"))) {
826 snprintf(net->atlas_space,64*sizeof(char),"%s",sbuf);
827 } else {
828 snprintf(net->atlas_space,64*sizeof(char),"UNKNOWN");
829 }
830 break ;
831 case NI_ELEMENT_TYPE:
832 nel = (NI_element *)ngr->part[ip] ;
833 if (!strcmp(nel->name,"tract") || !strcmp(nel->name,"tracts")) {
834 if ((tt = NIel_2_Tracts(nel, &N_tracts))) {
835 char *be=NULL;
836 tbb = AppCreateBundle(tbb, N_tracts, tt);
837 tt = Free_Tracts(tt, N_tracts);
838 NI_GET_INT(nel,"Bundle_Tag",ei);
839 if (!NI_GOT) ei = -1;
840 if ((be = NI_get_attribute(nel,"Bundle_Ends"))) {
841 net = AppAddBundleToNetwork(net, &tbb, ei, -1, NULL,be);
842 } else {
843 // Sept 2014
844 snprintf( tb_ends, 128, "%03d<->%s", ei,"-1");
845 net = AppAddBundleToNetwork(net, &tbb, ei, -1, NULL,
846 tb_ends);
847 }
848 } else {
849 WARNING_message("Failed to interpret nel tract,"
850 " ignoring.\n");
851 }
852 } else {
853 WARNING_message("Don't know about nel %s\n", nel->name);
854 }
855 break;
856 default:
857 ERROR_message("Don't know what to make of this "
858 "group element, ignoring.");
859 break;
860 }
861 }
862 }
863 RETURN(net);
864 }
865
AppAddBundleToNetwork(TAYLOR_NETWORK * network,TAYLOR_BUNDLE ** tb,int tag,int alt_tag,THD_3dim_dataset * grid,char * EleName)866 TAYLOR_NETWORK *AppAddBundleToNetwork(TAYLOR_NETWORK *network,
867 TAYLOR_BUNDLE **tb, int tag, int alt_tag,
868 THD_3dim_dataset *grid, char *EleName)
869 {
870 TAYLOR_NETWORK *net=NULL;
871
872 ENTRY("AppAddBundleToNetwork");
873
874 if (!tb) RETURN(net);
875
876 if (!network) {
877 net = (TAYLOR_NETWORK *)calloc(1,sizeof(TAYLOR_NETWORK));
878 net->N_allocated = -1;
879 net->N_points_private = -1;
880 if (grid) {
881 snprintf(net->atlas_space,64*sizeof(char),"%s", grid->atlas_space);
882 } else {
883 snprintf(net->atlas_space,64*sizeof(char),"UNKNOWN");
884 }
885 } else {
886 net = network;
887 net->N_points_private = -1; /* uninitialize so that it gets reset later */
888 }
889
890 if (net->N_allocated <= 0 || net->N_tbv >= net->N_allocated) {
891 net->N_allocated += 100;
892 net->tbv = (TAYLOR_BUNDLE **)realloc( net->tbv,
893 net->N_allocated*sizeof(TAYLOR_BUNDLE *));
894 net->bundle_tags = (int *)realloc(net->bundle_tags,
895 sizeof(int)*net->N_allocated);
896 net->bundle_alt_tags = (int *)realloc(net->bundle_alt_tags,
897 sizeof(int)*net->N_allocated);
898 }
899
900 /* PT: Let's chat about the location of bundle_ends and tags soon */
901 if (EleName) (*tb)->bundle_ends = strdup(EleName);
902 net->tbv[net->N_tbv] = *tb; *tb = NULL;
903 net->bundle_tags[net->N_tbv] = tag;
904 net->bundle_alt_tags[net->N_tbv] = alt_tag;
905 ++net->N_tbv;
906
907 RETURN(net);
908 }
909
Write_NI_Network(NI_group * ngr,char * name,char * mode)910 int Write_NI_Network(NI_group *ngr, char *name, char *mode)
911 {
912 char *nameout=NULL;
913 NI_stream ns;
914
915 ENTRY("Write_NI_Network");
916
917 if (!mode) mode = "NI_fast_binary";
918
919 /* be sure to init for tract datum */
920 if (get_NI_tract_type() < 0) {
921 ERROR_message("Misere!");
922 RETURN(0);
923 }
924 if (!name) name = "no_name";
925
926 nameout = (char *)calloc(strlen(name)+35, sizeof(char));
927 strcpy(nameout, "file:");
928 strcat(nameout,name);
929 nameout = without_afni_filename_extension(nameout);
930 strcat(nameout,".niml.tract");
931
932 ns = NI_stream_open(nameout, "w");
933 if (!ns) {
934 ERROR_message("Failed to open NI stream %s for writing.", nameout);
935 RETURN(0);
936 }
937
938 if (tract_verb) {
939 fprintf(stderr,"About to write %s in mode %s...", nameout, mode);
940 }
941 if (strcasestr(mode,"text")) {
942 NI_write_element( ns , ngr , NI_TEXT_MODE ) ;
943 } else {
944 NI_write_element( ns , ngr , NI_BINARY_MODE ) ;
945 }
946 if (tract_verb) {
947 fprintf(stderr," Done.\n");
948 }
949 NI_stream_close(ns); ns=NULL;
950 free(nameout);
951
952 RETURN(1);
953 }
954
Write_Bundle(TAYLOR_BUNDLE * tb,char * name,char * mode)955 int Write_Bundle(TAYLOR_BUNDLE *tb, char *name, char *mode)
956 {
957 TAYLOR_NETWORK *net=NULL;
958 int rval=0;
959
960 ENTRY("Write_Bundle");
961
962 if (!name) name = "no_name_jack";
963 if (!tb) RETURN(0);
964
965 net = (TAYLOR_NETWORK *)calloc(1,sizeof(TAYLOR_NETWORK));
966 net->tbv = (TAYLOR_BUNDLE**)calloc(1,sizeof(TAYLOR_BUNDLE*));
967 net->bundle_tags = (int *)calloc(1,sizeof(int));
968 net->bundle_alt_tags = (int *)calloc(1,sizeof(int));
969 net->tbv[0]=tb;
970 net->bundle_tags[0]=-1;
971 net->bundle_alt_tags[0]=-1;
972 net->N_tbv=1;
973
974 rval = Write_Network(net, name, mode);
975
976 net->tbv[0]=0; net->N_tbv=0;
977 Free_Network(net);
978 RETURN(rval);
979 }
980
Write_Network(TAYLOR_NETWORK * net,char * name,char * mode)981 int Write_Network(TAYLOR_NETWORK *net,
982 char *name, char *mode)
983 {
984 NI_group *ngr=NULL;
985 int rval=0;
986
987 ENTRY("Write_Network");
988 if (!name) name = "no_name_jack";
989 if (!net) RETURN(0);
990
991 if (!mode) mode = "NI_fast";
992 if (net->N_tbv > 1 && !strcasestr(mode,"NI_fast")) {
993 ERROR_message("Cannot write more than one bundle in slow mode");
994 RETURN(0);
995 }
996 if (strcasestr(mode,"NI_fast")) {
997 ngr = Network_2_NIgr(net, 1);
998 } else if (strcasestr(mode,"NI_slow")) {
999 ngr = Network_2_NIgr(net, 0);
1000 } else {
1001 ERROR_message("Stop making bad choices! %s\n",mode);
1002 RETURN(0);
1003 }
1004
1005 rval = Write_NI_Network(ngr, name, mode);
1006 NI_free_element(ngr); ngr=NULL;
1007
1008 RETURN(rval);
1009 }
1010
Read_NI_Network(char * name)1011 NI_group * Read_NI_Network(char *name)
1012 {
1013
1014 NI_stream ns;
1015 NI_group *ngr=NULL;
1016 char *nameout=NULL;
1017
1018 ENTRY("Read_NI_Network");
1019
1020 /* be sure to init for tract datum */
1021 if (get_NI_tract_type() < 0) {
1022 ERROR_message("Misere!");
1023 RETURN(ngr);
1024 }
1025
1026 if (!name) RETURN(ngr);
1027
1028 if (strcmp(name,"file:")) { /* have regular file name */
1029 if (THD_is_file(name)) {
1030 nameout = (char *)calloc(strlen(name)+35, sizeof(char));
1031 sprintf(nameout, "file:%s",name);
1032 } else {
1033 nameout = (char *)calloc(strlen(name)+35, sizeof(char));
1034 name = without_afni_filename_extension(name);
1035 sprintf(nameout,"%s.niml.tract", name);
1036 if (THD_is_file(nameout)) {
1037 sprintf(nameout, "file:%s.niml.tract", name);
1038 } else {
1039 ERROR_message("Cannot find %s\n", name);
1040 RETURN(ngr);
1041 }
1042 }
1043 }
1044
1045 ns = NI_stream_open(nameout,"r");
1046 if (!ns) RETURN(ngr);
1047 if (get_tract_verb()) fprintf(stderr,"About to read %s ...", nameout);
1048 ngr = NI_read_element( ns, 1) ;
1049 if (get_tract_verb()) fprintf(stderr," Done.\n");
1050
1051 NI_stream_close(ns);
1052 RETURN(ngr);
1053 }
1054
Read_Network(char * name)1055 TAYLOR_NETWORK * Read_Network(char *name)
1056 {
1057 NI_group *ngr=NULL;
1058 TAYLOR_NETWORK *net=NULL;
1059
1060 ENTRY("Read_Network");
1061
1062 if (!name) RETURN(net);
1063
1064 if (!(ngr = Read_NI_Network(name))) {
1065 ERROR_message("Failed to read NI_Bundle %s\n", name);
1066 RETURN(net);
1067 }
1068
1069 if (!(net = NIgr_2_Network(ngr))) {
1070 ERROR_message("Failed to turn group element to bundle %s\n", name);
1071 NI_free_element(ngr); ngr = NULL;
1072 RETURN(net);
1073 }
1074
1075 NI_free_element(ngr); ngr = NULL;
1076
1077 RETURN(net);
1078 }
1079
NI_getTractAlgOpts_M(NI_element * nel,float * MinFA,float * MaxAngDeg,float * MinL,int * SeedPerV)1080 int NI_getTractAlgOpts_M(NI_element *nel, float *MinFA, float *MaxAngDeg,
1081 float *MinL, int *SeedPerV)
1082 {
1083 char *atr=NULL;
1084
1085 ENTRY("NI_getTractAlgOpts");
1086 if (!nel) RETURN(1);
1087
1088 if (MinFA && (atr=NI_get_attribute(nel,"Thresh_FA"))) {
1089 *MinFA = (float)strtod(atr,NULL);
1090 }
1091 if (MaxAngDeg && (atr=NI_get_attribute(nel,"Thresh_ANG"))) {
1092 *MaxAngDeg = (float)strtod(atr,NULL);
1093 }
1094 if (MinL && (atr=NI_get_attribute(nel,"Thresh_Len"))) {
1095 *MinL = (float)strtod(atr,NULL);
1096 }
1097 if (SeedPerV && (atr=NI_get_attribute(nel,"Nseed_X"))) {
1098 SeedPerV[0] = (int)strtod(atr,NULL);
1099 }
1100 if (SeedPerV && (atr=NI_get_attribute(nel,"Nseed_Y"))) {
1101 SeedPerV[1] = (int)strtod(atr,NULL);
1102 }
1103 if (SeedPerV && (atr=NI_get_attribute(nel,"Nseed_Z"))) {
1104 SeedPerV[2] = (int)strtod(atr,NULL);
1105 }
1106
1107 RETURN(0);
1108 }
1109
NI_setTractAlgOpts_M(NI_element * nel,float * MinFA,float * MaxAngDeg,float * MinL,int * SeedPerV)1110 NI_element * NI_setTractAlgOpts_M(NI_element *nel, float *MinFA,
1111 float *MaxAngDeg, float *MinL,
1112 int *SeedPerV)
1113 {
1114 ENTRY("NI_setTractAlgOpts");
1115
1116 if (!nel) nel = NI_new_data_element ("TRACK_opts",0);
1117
1118 if (MinFA ) {
1119 NI_SETA_FLOAT(nel,"Thresh_FA",*MinFA);
1120 }
1121 if (MaxAngDeg) {
1122 NI_SETA_FLOAT(nel,"Thresh_ANG",*MaxAngDeg);
1123 }
1124 if (MinL) {
1125 NI_SETA_FLOAT(nel,"Thresh_Len",*MinL);
1126 }
1127 if (SeedPerV) {
1128 NI_SETA_INT(nel,"Nseed_X",SeedPerV[0]);
1129 NI_SETA_INT(nel,"Nseed_Y",SeedPerV[1]);
1130 NI_SETA_INT(nel,"Nseed_Z",SeedPerV[2]);
1131 }
1132
1133 RETURN(nel);
1134 }
1135
1136
ReadTractAlgOpts_M(char * fname)1137 NI_element * ReadTractAlgOpts_M(char *fname)
1138 {
1139 NI_stream ns=NULL;
1140 NI_element *nel=NULL;
1141 float MinFA, MaxAngDeg, MinL;
1142 int SeedPerV[3];
1143 char *strm=NULL;
1144 FILE *fin4=NULL;
1145
1146 ENTRY("ReadTractAlgOpts");
1147
1148 if (!fname || !THD_is_file(fname)) RETURN(NULL);
1149
1150 if (STRING_HAS_SUFFIX(fname,".niml.opts")) {
1151 strm = (char *)calloc(strlen(fname)+20, sizeof(char));
1152 sprintf(strm,"file:%s",fname);
1153 if (!(ns = NI_stream_open( strm , "r" ))) {
1154 ERROR_message("Failed to open %s\n", strm);
1155 free(strm); RETURN(NULL);
1156 }
1157 if (!(nel = NI_read_element( ns , 2 ))) {
1158 ERROR_message("Failed to read element from \n", strm);
1159 free(strm); RETURN(NULL);
1160 }
1161 NI_stream_close(ns); free(strm); strm = NULL;
1162 } else {
1163 // Opening/Reading in FACT params
1164 if( (fin4 = fopen(fname, "r")) == NULL) {
1165 fprintf(stderr, "Error opening file %s.",fname);
1166 RETURN(NULL);
1167 }
1168
1169 if( !(fscanf(fin4, "%f %f %f %d %d %d",
1170 &MinFA,&MaxAngDeg,&MinL,&SeedPerV[0],&SeedPerV[1],
1171 &SeedPerV[2]))){
1172 fprintf(stderr, "Error reading parameter files.");
1173 RETURN(NULL);
1174 }
1175
1176 fclose(fin4);
1177 if (!(nel =
1178 NI_setTractAlgOpts_M(NULL, &MinFA, &MaxAngDeg, &MinL,
1179 SeedPerV))){
1180 ERROR_message("Failed to get options");
1181 RETURN(NULL);
1182 }
1183 }
1184
1185 RETURN(nel);
1186 }
1187
1188 // THIS one doesn't need to change with new format for MULTI/HARDI update
WriteTractAlgOpts(char * fname,NI_element * nel)1189 int WriteTractAlgOpts(char *fname, NI_element *nel)
1190 {
1191 char *strm=NULL;
1192 NI_stream ns=NULL;
1193
1194 ENTRY("WriteTractAlgOpts");
1195
1196 if (!nel) {
1197 fprintf(stderr, "NULL nel\n");
1198 RETURN(1);
1199 }
1200 if (fname) {
1201 strm = (char *)calloc(strlen(fname)+20, sizeof(char));
1202 if (STRING_HAS_SUFFIX(fname,".niml.opts")) {
1203 sprintf(strm,"file:%s",fname);
1204 } else {
1205 sprintf(strm,"file:%s.niml.opts",fname);
1206 }
1207 } else {
1208 strm = (char *)calloc(20, sizeof(char));
1209 sprintf(strm,"fd:1");
1210 }
1211 if (!(ns = NI_stream_open( strm , "w" ))) {
1212 ERROR_message("Failed to open %s\n", strm);
1213 free(strm); RETURN(1);
1214 }
1215 NI_write_element(ns,nel,NI_TEXT_MODE);
1216 NI_stream_close(ns); free(strm); strm = NULL;
1217 RETURN(0);
1218 }
1219
NI_getProbTractAlgOpts_M(NI_element * nel,float * MinFA,float * MaxAngDeg,float * MinL,float * NmNsFr,int * Nseed,int * Nmonte)1220 int NI_getProbTractAlgOpts_M(NI_element *nel, float *MinFA, float *MaxAngDeg,
1221 float *MinL, float *NmNsFr, int *Nseed,
1222 int *Nmonte)
1223 {
1224 char *atr=NULL;
1225
1226 ENTRY("NI_getProbTractAlgOpts");
1227 if (!nel) RETURN(1);
1228
1229 if (MinFA && ( (atr=NI_get_attribute(nel,"Thresh_FA")) ||
1230 (atr=NI_get_attribute(nel,"MinFA"))) ) {
1231 *MinFA = (float)strtod(atr,NULL);
1232 }
1233 if (MaxAngDeg && ( (atr=NI_get_attribute(nel,"Thresh_ANG")) ||
1234 (atr=NI_get_attribute(nel,"MaxAng"))) ) {
1235 *MaxAngDeg = (float)strtod(atr,NULL);
1236 }
1237 if (MinL && ( (atr=NI_get_attribute(nel,"Thresh_Len")) ||
1238 (atr=NI_get_attribute(nel,"MinL"))) ) {
1239 *MinL = (float)strtod(atr,NULL);
1240 }
1241
1242 if (NmNsFr && ( (atr=NI_get_attribute(nel,"Thresh_Frac")) ||
1243 (atr=NI_get_attribute(nel,"MinHitFr"))) ) {
1244 *NmNsFr = (float)strtod(atr,NULL);
1245 }
1246 if (Nseed && ( (atr=NI_get_attribute(nel,"Nseed_Vox")) ||
1247 (atr=NI_get_attribute(nel,"Nseed"))) ) {
1248 *Nseed = (int)strtod(atr,NULL);
1249 }
1250 if (Nmonte && (atr=NI_get_attribute(nel,"Nmonte"))) {
1251 *Nmonte = (int)strtod(atr,NULL);
1252 }
1253 RETURN(0);
1254 }
1255
NI_setProbTractAlgOpts_M(NI_element * nel,float * MinFA,float * MaxAngDeg,float * MinL,float * NmNsFr,int * Nseed,int * Nmonte)1256 NI_element * NI_setProbTractAlgOpts_M(NI_element *nel, float *MinFA,
1257 float *MaxAngDeg, float *MinL,
1258 float *NmNsFr, int *Nseed,
1259 int *Nmonte)
1260 {
1261 ENTRY("NI_setProbTractAlgOpts_M");
1262
1263 if (!nel) nel = NI_new_data_element ("PROBTRACK_opts",0);
1264
1265 if (MinFA ) {
1266 NI_SETA_FLOAT(nel,"Thresh_FA",*MinFA);
1267 }
1268 if (MaxAngDeg) {
1269 NI_SETA_FLOAT(nel,"Thresh_ANG",*MaxAngDeg);
1270 }
1271 if (MinL) {
1272 NI_SETA_FLOAT(nel,"Thresh_Len",*MinL);
1273 }
1274 if (NmNsFr) {
1275 NI_SETA_FLOAT(nel,"Thresh_Frac",*NmNsFr);
1276 }
1277 if (Nseed) {
1278 NI_SETA_INT(nel,"Nseed_Vox",*Nseed);
1279 }
1280 if (Nmonte) {
1281 NI_SETA_INT(nel,"Nmonte",*Nmonte);
1282 }
1283
1284 RETURN(nel);
1285 }
1286
ReadProbTractAlgOpts_M(char * fname)1287 NI_element * ReadProbTractAlgOpts_M(char *fname)
1288 {
1289 NI_stream ns=NULL;
1290 NI_element *nel=NULL;
1291 float MinFA, MaxAngDeg, MinL,NmNsFr;
1292 int Nseed, Nmonte;
1293 char *strm=NULL;
1294 FILE *fin4=NULL;
1295
1296 ENTRY("ReadProbTractAlgOpts");
1297
1298 if (!fname || !THD_is_file(fname)) RETURN(NULL);
1299
1300 if (STRING_HAS_SUFFIX(fname,".niml.opts")) {
1301 strm = (char *)calloc(strlen(fname)+20, sizeof(char));
1302 sprintf(strm,"file:%s",fname);
1303 if (!(ns = NI_stream_open( strm , "r" ))) {
1304 ERROR_message("Failed to open %s\n", strm);
1305 free(strm); RETURN(NULL);
1306 }
1307 if (!(nel = NI_read_element( ns , 2 ))) {
1308 ERROR_message("Failed to read element from \n", strm);
1309 free(strm); RETURN(NULL);
1310 }
1311 NI_stream_close(ns); free(strm); strm = NULL;
1312 } else {
1313 // Opening/Reading in FACT params
1314 if( (fin4 = fopen(fname, "r")) == NULL) {
1315 fprintf(stderr, "Error opening file %s.",fname);
1316 RETURN(NULL);
1317 }
1318
1319 if( !(fscanf(fin4, "%f %f %f %f %d %d",
1320 &MinFA,&MaxAngDeg,&MinL,&NmNsFr,&Nseed,&Nmonte))){
1321 fprintf(stderr, "Error reading parameter files.");
1322 RETURN(NULL);
1323 }
1324 fclose(fin4);
1325 //printf("%f %f %f %f %d %d",
1326 // MinFA,MaxAngDeg,MinL,NmNsFr,Nseed,Nmonte);
1327
1328 if (!(nel =
1329 NI_setProbTractAlgOpts_M(NULL, &MinFA, &MaxAngDeg, &MinL,
1330 &NmNsFr,&Nseed,&Nmonte))){
1331 ERROR_message("Failed to get options");
1332 RETURN(NULL);
1333 }
1334 }
1335
1336 RETURN(nel);
1337 }
1338
SimpleWriteDetNetTr_M(int N_HAR,FILE * file,int *** idx,THD_3dim_dataset ** PARS,int PAR_BOT,int PAR_TOP,float ** loc,int ** locI,int len,int * TV,int * Dim,float * Ledge)1339 int SimpleWriteDetNetTr_M(int N_HAR, FILE *file, int ***idx,
1340 THD_3dim_dataset **PARS,
1341 int PAR_BOT, int PAR_TOP,
1342 float **loc, int **locI, int len,
1343 int *TV, int *Dim, float *Ledge)
1344 { // for trackvis
1345 int m,aa,bb;
1346 int READS_in;
1347 float READS_fl;
1348
1349 ENTRY("SimpleWriteDetNetTr");
1350
1351 // first write len of tr
1352 READS_in = len;
1353 fwrite(&READS_in,sizeof(READS_in),1,file);
1354
1355 // then track, coors and props
1356 for( m=0 ; m<len ; m++ ) {
1357 // recenter phys loc for trackvis, if nec...
1358 for( aa=0 ; aa<3 ; aa++ ) {
1359 READS_fl = loc[m][aa];
1360 if(!TV[aa])
1361 READS_fl = Ledge[aa]*Dim[aa]-READS_fl;
1362 fwrite(&READS_fl,sizeof(READS_fl),1,file);
1363 }
1364 bb=idx[locI[m][0]][locI[m][1]][locI[m][2]];
1365 if(N_HAR) { // hardi, single param
1366 READS_fl = THD_get_voxel(PARS[PAR_BOT], bb, 0);
1367 fwrite(&READS_fl,sizeof(READS_fl),1,file);
1368 }
1369 else
1370 for( aa=1 ; aa<4 ; aa++ ) {
1371 READS_fl = THD_get_voxel(PARS[aa], bb, 0);
1372 fwrite(&READS_fl,sizeof(READS_fl),1,file);
1373 }
1374 }
1375
1376 RETURN(1);
1377 }
1378
Free_Insta_Tract_Setup(INSTA_TRACT_SETUP * ITS)1379 int Free_Insta_Tract_Setup(INSTA_TRACT_SETUP *ITS)
1380 {
1381 ENTRY("Free_Insta_Tract_Setup");
1382
1383 if (!ITS) RETURN(0);
1384
1385 if (ITS->grid) DSET_delete(ITS->grid);
1386 ITS->grid = NULL;
1387
1388 /* Do not delte ITS , leave it to calling function */
1389
1390 RETURN(1);
1391 }
1392
1393 /* Create brandnew, or wipe clean existing ITS */
New_Insta_Tract_Setup(INSTA_TRACT_SETUP * ITS)1394 INSTA_TRACT_SETUP *New_Insta_Tract_Setup(INSTA_TRACT_SETUP *ITS)
1395 {
1396 ENTRY("New_Insta_Tract_Setup");
1397
1398 if (!ITS) ITS = (INSTA_TRACT_SETUP *)calloc(1,sizeof(INSTA_TRACT_SETUP));
1399 else Free_Insta_Tract_Setup(ITS);
1400
1401 /* Put any initialization here ... */
1402
1403 RETURN(ITS);
1404 }
1405
1406
1407 // *****************-> NIMLly reading in DTI input <-*********************
1408
ReadDTI_inputs(char * fname)1409 NI_element * ReadDTI_inputs(char *fname)
1410 {
1411 NI_stream ns=NULL;
1412 NI_element *nel=NULL;
1413 char *strm=NULL;
1414 FILE *fin4=NULL;
1415
1416 ENTRY("ReadDTI_inputs");
1417
1418 if (!fname || !THD_is_file(fname)) RETURN(NULL);
1419
1420 if (STRING_HAS_SUFFIX(fname,".niml.opts")) {
1421 strm = (char *)calloc(strlen(fname)+20, sizeof(char));
1422 sprintf(strm,"file:%s",fname);
1423 if (!(ns = NI_stream_open( strm , "r" ))) {
1424 ERROR_message("Failed to open %s\n", strm);
1425 free(strm); RETURN(NULL);
1426 }
1427 if (!(nel = NI_read_element( ns , 2 ))) {
1428 ERROR_message("Failed to read element from \n", strm);
1429 free(strm); RETURN(NULL);
1430 }
1431 NI_stream_close(ns); free(strm); strm = NULL;
1432 } else {
1433 ERROR_message("Failed to get DTI inputs from %s",fname);
1434 RETURN(NULL);
1435 }
1436
1437 RETURN(nel);
1438 }
1439
1440
NI_getDTI_inputs(NI_element * nel,char ** NameVECT,char * NameXF,char ** NameSCAL,char ** NamePLUS,int * extrafile,int * pars_top)1441 int NI_getDTI_inputs( NI_element *nel,
1442 char **NameVECT,
1443 char *NameXF,
1444 char **NameSCAL,
1445 char **NamePLUS,
1446 int *extrafile, int *pars_top)
1447 {
1448 char *atr="NONAME";
1449 char tmp[THD_MAX_PREFIX];
1450 int i;
1451 int ct_scal = 1; // 1 -> start with space for 'extrafile'
1452
1453 ENTRY("NI_getDTI_inputs");
1454 if (!nel) RETURN(1);
1455
1456 atr = (char *)calloc(100, sizeof(char));
1457 if( atr == NULL ) {
1458 fprintf(stderr, "\n\n MemAlloc failure.\n\n");
1459 exit(126);
1460 }
1461
1462 // for vectors
1463 for( i=0 ; i<N_DTI_VECT ; i++ ) {
1464 sprintf(tmp, "dti_%s", DTI_VECT_LABS[i]);
1465 if (NameVECT[i] && (atr=NI_get_attribute(nel,tmp))) {
1466 snprintf(NameVECT[i], N_CHAR_PATH, "%s", atr);
1467 }
1468 }
1469
1470 // for scalars
1471 for( i=0 ; i<N_DTI_SCAL ; i++ ) {
1472 sprintf(tmp, "dti_%s", DTI_SCAL_LABS[i]);
1473 if (NameSCAL[i] && (atr=NI_get_attribute(nel,tmp))) {
1474 snprintf(NameSCAL[i], N_CHAR_PATH, "%s", atr);
1475 ct_scal++;
1476 }
1477 }
1478
1479 // for extra scalar
1480 sprintf(tmp, "dti_%s", DTI_XTRA_LABS[0]);
1481 if (NameXF && (atr=NI_get_attribute(nel,tmp))) {
1482 snprintf(NameXF, N_CHAR_PATH, "%s", atr);
1483 *extrafile = 1;
1484 }
1485 else
1486 NameXF = NULL;
1487
1488
1489 // allow up to four extra files
1490 for( i=0 ; i<N_DTI_PLUS ; i++ ) {
1491 sprintf(tmp, "dti_%s", DTI_PLUS_LABS[i]);
1492 if (NamePLUS[i] && (atr=NI_get_attribute(nel,tmp))) {
1493 snprintf(NamePLUS[i], N_CHAR_PATH, "%s", atr);
1494 ct_scal++;
1495 }
1496 else
1497 snprintf(NamePLUS[i], N_CHAR_PATH, "%s", "\0");
1498 }
1499
1500 *pars_top = ct_scal; // 2 + 3 + ct_ex; // RD and extra; FA,MD;extras
1501 INFO_message(" ct_scal: %d atr:%s ", ct_scal, atr);
1502
1503 RETURN(0);
1504 }
1505
1506
SUMA_Taylor_Bundle_Info(TAYLOR_BUNDLE * tb,int show_maxu)1507 char *SUMA_Taylor_Bundle_Info(TAYLOR_BUNDLE *tb, int show_maxu)
1508 {
1509 static char FuncName[]={"SUMA_Taylor_Bundle_Info"};
1510 int show_max;
1511 int ii=0;
1512 char stmp[64];
1513 char *s=NULL;
1514 SUMA_STRING *SS=NULL;
1515
1516 SUMA_ENTRY;
1517
1518 SS = SUMA_StringAppend(NULL, NULL);
1519
1520 if (!tb) {
1521 SUMA_StringAppend(SS,"NULL bundle pointer");
1522 } else {
1523 SUMA_StringAppend_va(SS, "Bundle has %d tracts\n", tb->N_tracts);
1524 if ((show_maxu < 0) || (tb->N_tracts < show_maxu)) show_max = tb->N_tracts;
1525 else if (show_maxu == 0) show_max = (tb->N_tracts < 5) ? tb->N_tracts : 5;
1526 else show_max = show_maxu < tb->N_tracts ? show_maxu : tb->N_tracts;
1527
1528 s = NULL;
1529 for (ii=0; ii<show_max; ++ii) {
1530 snprintf(stmp, 62," Bun.Trc %d ++> ", ii);
1531 s = SUMA_append_replace_string(s,
1532 SUMA_Taylor_Tract_Info(tb->tracts+ii, show_maxu),stmp,2);
1533 }
1534 SUMA_StringAppend_va(SS,s); SUMA_ifree(s);
1535 if (show_max < tb->N_tracts) {
1536 int pl = (tb->N_tracts-show_max > 1);
1537 SUMA_StringAppend_va(SS,
1538 " ... %d tract%sremain%s in bundle.\n",
1539 tb->N_tracts-show_max, pl ? "s ":" ", pl ? "":"s");
1540 }
1541 }
1542
1543 SUMA_SS2S(SS, s);
1544 SUMA_RETURN(s);
1545 }
1546
1547
1548
SUMA_Taylor_Tract_Info(TAYLOR_TRACT * tt,int show_maxu)1549 char *SUMA_Taylor_Tract_Info(TAYLOR_TRACT *tt, int show_maxu)
1550 {
1551 static char FuncName[]={"SUMA_Taylor_Tract_Info"};
1552 int show_max;
1553 int ii=0;
1554 char *s=NULL;
1555 SUMA_STRING *SS=NULL;
1556
1557 SUMA_ENTRY;
1558
1559 SS = SUMA_StringAppend(NULL, NULL);
1560
1561 if (!tt) {
1562 SUMA_StringAppend(SS,"NULL tract pointer");
1563 } else {
1564 SUMA_StringAppend_va(SS, " track id %d, Npts=%d\n",
1565 tt->id, TRACT_NPTS(tt));
1566 if (show_maxu < 0) show_max = TRACT_NPTS(tt);
1567 else if (show_maxu == 0)
1568 show_max = (TRACT_NPTS(tt) < 5) ? TRACT_NPTS(tt) : 5;
1569 else show_max = show_maxu < TRACT_NPTS(tt) ? show_maxu : TRACT_NPTS(tt);
1570 for (ii=0; ii<show_max; ++ii) {
1571 SUMA_StringAppend_va(SS, " %d %f %f %f\n",
1572 ii, tt->pts[3*ii], tt->pts[3*ii+1],tt->pts[3*ii+2]);
1573 }
1574 if (show_max < TRACT_NPTS(tt)) {
1575 int pl = (TRACT_NPTS(tt)-show_max > 1);
1576 SUMA_StringAppend_va(SS,
1577 " ... %d point%sremain%s in tract.\n",
1578 TRACT_NPTS(tt)-show_max, pl ? "s ":" ", pl ? "":"s");
1579 }
1580 }
1581
1582 SUMA_SS2S(SS, s);
1583 SUMA_RETURN(s);
1584 }
1585
1586
1587 /* Both of SUMA_Network_N_points and SUMA_Network_N_tracts
1588 are copies of the one in ptaylor/TrackIO.c
1589 The other ones should be used as soon as
1590 we resolve afni's dependency on ptaylor/ */
SUMA_Network_N_points(TAYLOR_NETWORK * net,byte recalc)1591 int SUMA_Network_N_points(TAYLOR_NETWORK *net, byte recalc)
1592 {
1593 int nb, ib=0, nn=-1, it;
1594
1595 if (!net) return(-1);
1596 if (!recalc && net->N_points_private > 0) return(net->N_points_private);
1597 nn = 0;
1598 for (ib=0; ib<net->N_tbv; ++ib) {
1599 if (net->tbv[ib]) {
1600 nb = 0;
1601 for (it=0; it<net->tbv[ib]->N_tracts; ++it) {
1602 nb += net->tbv[ib]->tracts[it].N_pts3;
1603 }
1604 net->tbv[ib]->N_points_private = nb/3;
1605 nn += nb;
1606 }
1607 }
1608 nn /= 3;
1609 net->N_points_private = nn;
1610 return(nn);
1611 }
1612
SUMA_Network_N_tracts(TAYLOR_NETWORK * net,byte recalc)1613 int SUMA_Network_N_tracts(TAYLOR_NETWORK *net, byte recalc)
1614 {
1615 int nb, ib=0, nn=-1, it;
1616
1617 if (!net) return(-1);
1618 if (!recalc && net->N_tracts_private > 0) return(net->N_tracts_private);
1619 nn = 0;
1620 for (ib=0; ib<net->N_tbv; ++ib) {
1621 if (net->tbv[ib]) {
1622 nn += net->tbv[ib]->N_tracts;
1623 }
1624 }
1625 net->N_tracts_private = nn;
1626 return(nn);
1627 }
1628
1629
SUMA_Taylor_Network_Info(TAYLOR_NETWORK * net,int show_maxu,int show_maxub)1630 char * SUMA_Taylor_Network_Info(TAYLOR_NETWORK *net,
1631 int show_maxu, int show_maxub)
1632 {
1633 static char FuncName[]={"SUMA_Taylor_Network_Info"};
1634 TAYLOR_BUNDLE *tb=NULL;
1635 int show_max;
1636 char stmp[64];
1637 int ii=0;
1638 char *s=NULL;
1639 SUMA_STRING *SS=NULL;
1640
1641 SUMA_ENTRY;
1642
1643 SS = SUMA_StringAppend(NULL, NULL);
1644
1645 if (!net) {
1646 SUMA_StringAppend(SS,"NULL network pointer");
1647 } else {
1648 SUMA_StringAppend_va(SS," Network has %d bundles, %d tracts, %d points\n",
1649 net->N_tbv, SUMA_Network_N_tracts(net,1),
1650 SUMA_Network_N_points(net, 1));
1651 if (show_maxu < 0) show_max = net->N_tbv;
1652 else if (show_maxu == 0) show_max = (net->N_tbv < 5) ? net->N_tbv : 5;
1653 else show_max = show_maxu < net->N_tbv ? show_maxu : net->N_tbv;
1654
1655 s = NULL;
1656 for (ii=0; ii<show_max; ++ii) {
1657 snprintf(stmp, 62, " Net.Bun. %d --> ", ii);
1658 s = SUMA_append_replace_string(s,
1659 SUMA_Taylor_Bundle_Info(net->tbv[ii], show_maxub),stmp,2);
1660 }
1661 SUMA_StringAppend_va(SS, s); SUMA_ifree(s);
1662 if (show_max < net->N_tbv) {
1663 int pl = (net->N_tbv-show_max > 1);
1664 SUMA_StringAppend_va(SS,
1665 "... %d bundle%sremain%s in network.\n",
1666 net->N_tbv-show_max, pl ? "s ":" ", pl ? "":"s");
1667 }
1668
1669 }
1670
1671 SUMA_SS2S(SS, s);
1672 SUMA_RETURN(s);
1673 }
1674 //OLD
1675 /*int NI_getDTI_inputs( NI_element *nel,
1676 char **NameVEC,
1677 char *NameXF,
1678 char **NameSCAL,
1679 char **NameP,
1680 int *extrafile, int *pars_top)
1681 {
1682 char *atr=NULL;
1683 int ct_ex = 0;
1684
1685 ENTRY("NI_getDTI_inputs");
1686 if (!nel) RETURN(1);
1687
1688 if (NameVEC[0] && (atr=NI_get_attribute(nel,"dti_V1"))) {
1689 snprintf(NameVEC[0],100,"%s", atr);
1690 }
1691 if (NameVEC[1] && (atr=NI_get_attribute(nel,"dti_V2"))) {
1692 snprintf(NameVEC[1],100,"%s", atr);
1693 }
1694 if (NameVEC[2] && (atr=NI_get_attribute(nel,"dti_V3"))) {
1695 snprintf(NameVEC[2],100,"%s", atr);
1696 }
1697 if (NameXF && (atr=NI_get_attribute(nel,"dti_XF"))) {
1698 snprintf(NameXF,100,"%s", atr); *extrafile = 1;
1699 }
1700 else
1701 NameXF = NULL;
1702 if (NameSCAL[0] && (atr=NI_get_attribute(nel,"dti_FA"))) {
1703 snprintf(NameSCAL[0],100,"%s", atr);
1704 }
1705 if (NameSCAL[1] && (atr=NI_get_attribute(nel,"dti_MD"))) {
1706 snprintf(NameSCAL[1],100,"%s", atr);
1707 }
1708 if (NameSCAL[2] && (atr=NI_get_attribute(nel,"dti_L1"))) {
1709 snprintf(NameSCAL[2],100,"%s", atr);
1710 }
1711 // allow up to four extra files
1712 if (NameP[0] && (atr=NI_get_attribute(nel,"dti_P1"))) {
1713 snprintf(NameP[0],100,"%s", atr); ct_ex++;
1714 }
1715 else snprintf(NameP[0],100,"%s", "\0");
1716 if (NameP[1] && (atr=NI_get_attribute(nel,"dti_P2"))) {
1717 snprintf(NameP[1],100,"%s", atr); ct_ex++;
1718 }
1719 else snprintf(NameP[1],100,"%s", "\0");
1720 if (NameP[2] && (atr=NI_get_attribute(nel,"dti_P3"))) {
1721 snprintf(NameP[2],100,"%s", atr); ct_ex++;
1722 }
1723 else snprintf(NameP[2],100,"%s", "\0");
1724 if (NameP[3] && (atr=NI_get_attribute(nel,"dti_P4"))) {
1725 snprintf(NameP[3],100,"%s", atr); ct_ex++;
1726 //printf("\nHERE! %d\n",ct_ex);
1727 }
1728 else snprintf(NameP[3],100,"%s", "\0");
1729
1730 *pars_top = 2 + 3 + ct_ex; // RD and extra; FA,MD;extras
1731
1732 //printf("\nLISTED: PARS_TOP=%d with ct_ex=%d\n", *pars_top, ct_ex);
1733 // printf("\nLISTED: NameSCAL[0]=%s\n", NameSCAL[0]);
1734
1735
1736 RETURN(0);
1737 }
1738
1739 */
1740