1 /*****************************************************************************
2    Major portions of this software are copyrighted by the Medical College
3    of Wisconsin, 1994-2000, and are released under the Gnu General Public
4    License, Version 2.  See the file README.Copyright for details.
5 ******************************************************************************/
6 
7 #include "mrilib.h"
8 #include "thd.h"
9 
10 /*---------------------------------------------------------------------
11    Given a directory name, read in all the datasets and make
12    a session data structure from them.
13 
14    [28 Jul 2003] Modified to use new THD_session struct, wherein all
15                  datasets are in one array (no anat/func distinction).
16 -----------------------------------------------------------------------*/
17 
THD_init_session(char * sessname)18 THD_session * THD_init_session( char *sessname )
19 {
20    THD_session            *sess ;
21    RwcPointer_array        *dblk_arrarr ;
22    THD_datablock_array    *dblk_arr ;
23    THD_3dim_dataset       *dset=NULL , *temp_dset;
24    THD_3dim_dataset_array *dset_arr ;
25 
26    int ibar , idset , iview  , nds , qview=-666 ;
27 
28 ENTRY("THD_init_session") ;
29 
30    /*-- sanity check --*/
31 
32    if( sessname == NULL || strlen(sessname) == 0 || !THD_is_directory(sessname) )
33      RETURN( NULL ) ;
34 
35    /*-- initialize session --*/
36 
37    sess         = myRwcNew( THD_session ) ;
38    sess->type   = SESSION_TYPE ;
39    sess->parent = NULL ;
40 #ifdef DEBUG_SESSIONS
41 fprintf(stderr, "blanking session\n");
42 #endif
43    BLANK_SESSION(sess) ;  /* null out all entries */
44 
45    /* save directory name, with a trailing slash */
46 
47    MCW_strncpy( sess->sessname , sessname , THD_MAX_NAME ) ;
48    iview = strlen(sess->sessname) ;
49    if( sess->sessname[iview-1] != '/' ){  /* tack trailing / onto sessname */
50      sess->sessname[iview]   = '/' ;
51      sess->sessname[iview+1] = '\0' ;
52    } else {
53      iview-- ;  /* iview now points to last non-NUL character in string */
54    }
55 
56    /* save last name from sessname */
57 #if 1
58    { char *env = my_getenv( "AFNI_SESSTRAIL" ) ; int tt = 0 ;
59      if( env != NULL ) tt = strtol(env,NULL,10) ;
60      env = THD_trailname(sess->sessname,tt) ;
61      tt = 1+strlen(env) - THD_MAX_NAME ; if( tt < 0 ) tt = 0 ;
62      strcpy( sess->lastname , env+tt ) ;
63    }
64 #else
65      for( iview-- ; iview >= 0 ; iview-- )
66        if( sess->sessname[iview] == '/' ) break ;
67      MCW_strncpy( sess->lastname, &(sess->sessname[iview+1]), THD_MAX_NAME ) ;
68 #endif
69 
70    /*-- override dataset 'view'??  [30 May 2013] --*/
71 
72    { char *env = my_getenv("AFNI_OVERRIDE_VIEW") ;
73      if( env != NULL ){
74             if( toupper(*env) == 'T' ) qview = VIEW_TALAIRACH_TYPE ;
75        else if( toupper(*env) == 'O' ) qview = VIEW_ORIGINAL_TYPE  ;
76        else
77          WARNING_message("AFNI_OVERRIDE_VIEW=%s is not understood",env) ;
78      }
79    }
80 
81    /*-- read all datablocks --*/
82 
83    dblk_arrarr = THD_init_alldir_datablocks( sess->sessname ) ;
84 
85    /*-- for each datablock array ... --*/
86 
87    for( ibar=0 ; ibar < dblk_arrarr->num ; ibar++ ){
88 
89       /*-- get the current array of datablocks --*/
90 
91       dblk_arr = (THD_datablock_array *) dblk_arrarr->ar[ibar] ;
92       if( dblk_arr == NULL || dblk_arr->num <= 0 ) continue ;    /* huh? */
93 
94       /*-- convert it into an array of datasets --*/
95 
96       dset_arr = THD_array_3dim_from_block( dblk_arr ) ;
97       if( dset_arr == NULL || dset_arr->num <= 0 ) continue ;
98 
99       /*-- place it into the next row of the THD_session [28 Jul 2003] --*/
100 
101       nds = sess->num_dsset ;
102 
103       if( nds >= THD_MAX_SESSION_SIZE ){   /* bad! */
104         fprintf(stderr,
105          "\n*** Session %s table overflow with dataset %s ***\n",
106              sessname , dset_arr->ar[0]->self_name) ;
107         for( idset=0 ; idset < dset_arr->num ; idset++ )
108           THD_delete_3dim_dataset( dset_arr->ar[idset] , False ) ;
109         FREE_3DARR(dset_arr) ;
110         continue ;  /* skip to next dblk_arr (ibar loop) */
111       }
112 
113       /*-- put each dataset into this row in its view place --*/
114 
115       for( idset=0 ; idset < dset_arr->num ; idset++ ){
116         dset  = dset_arr->ar[idset] ;
117         iview = dset->view_type ;
118         if( qview >= 0 ) iview = qview ;  /* 30 May 2013 */
119 
120         if( GET_SESSION_DSET(sess,nds,iview) != NULL ){  /* should never happen */
121           fprintf(stderr,
122            "\n*** Session %s has duplicate dataset views of %s ***\n",
123            sessname , dset->self_name) ;
124           THD_delete_3dim_dataset( dset , False ) ;
125         } else {
126 #ifdef DEBUG_SESSIONS
127 fprintf(stderr,"\nputting datasets into initial view \n");
128 #endif
129           SET_SESSION_DSET(dset, sess, nds, iview);  /* should always happen */
130         }
131       }
132 
133       sess->num_dsset ++ ;  /* one more row */
134 
135       FREE_3DARR(dset_arr) ;
136 
137    } /* end of loop over each datablock array (ibar) */
138 
139    /*-- throw away the datablock arrays at this point --*/
140 
141    STATUS("trashing dblk_arrarr") ;
142 
143    for( ibar=0 ; ibar < dblk_arrarr->num ; ibar++ ){
144      dblk_arr = (THD_datablock_array *) dblk_arrarr->ar[ibar] ;
145      FREE_DBARR( dblk_arr ) ;
146    }
147    FREE_XTARR( dblk_arrarr ) ;
148 
149    /*-- 06 Apr 2005: try to read NIfTI-1 files [KRH and RWC] --*/
150 
151    if( !AFNI_noenv("AFNI_NIFTI_DATASETS") ){
152      char *ename[2] , **fn_nifti ;
153      int num_nifti , ii ;
154 
155      STATUS("looking for NIFTI files") ;
156 
157      ename[0] = AFMALL(char, THD_MAX_NAME) ;
158      ename[1] = AFMALL(char, THD_MAX_NAME) ;
159      strcpy(ename[0],sess->sessname) ; strcat(ename[0],"*.nii") ;
160      strcpy(ename[1],sess->sessname) ; strcat(ename[1],"*.nii.gz") ;
161      MCW_file_expand( 2,ename , &num_nifti,&fn_nifti ) ;  /* find files */
162      free(ename[0]) ; free(ename[1]) ;
163 
164      if( num_nifti > 0 ){                               /* got some! */
165        STATUS("opening NIFTI files") ;
166        for( ii=0 ; ii < num_nifti ; ii++ ){             /* loop over files */
167          dset = THD_open_nifti( fn_nifti[ii] ) ;        /* try it on */
168          if( !ISVALID_DSET(dset) ) continue ;           /* doesn't fit? */
169          nds = sess->num_dsset ;
170          if( nds >= THD_MAX_SESSION_SIZE ){
171            fprintf(stderr,
172              "\n*** Session %s table overflow with NIfTI dataset %s ***\n",
173              sessname , fn_nifti[ii] ) ;
174            THD_delete_3dim_dataset( dset , False ) ;
175            break ; /* out of for loop */
176          }
177          iview = dset->view_type ;
178          SET_SESSION_DSET(dset, sess, nds, iview);
179          sess->num_dsset ++ ;
180        } /* end of loop over files */
181        MCW_free_expand( num_nifti , fn_nifti ) ;
182      } /* end of if we found NIFTI files */
183    }
184 
185    /*-- 27 Aug 2002: try to read any ANALYZE "datasets" here --*/
186 
187    if( !AFNI_noenv("AFNI_ANALYZE_DATASETS") ){
188      char *ename[2] , **fn_anlz ;
189      int num_anlz , ii , nee ;
190 #ifdef ALLOW_FSL_FEAT
191      int feat_exf=-1 , feat_hrs=-1 , feat_std=-1 ;
192      int feat_nds_start=sess->num_dsset ;
193 #endif
194 
195      STATUS("looking for ANALYZE files") ;
196 
197      ename[0] = AFMALL(char, THD_MAX_NAME) ;
198      strcpy(ename[0],sess->sessname) ; strcat(ename[0],"*.hdr") ;
199      nee = 1 ;
200 #ifdef ALLOW_FSL_FEAT
201      ename[1] = AFMALL(char, THD_MAX_NAME) ;
202      strcpy(ename[1],sess->sessname) ; strcat(ename[1],"stats/*stat*.hdr") ;
203      nee++ ;
204 #endif
205      MCW_file_expand( nee,ename , &num_anlz,&fn_anlz ) ;  /* find files */
206      for( ii=0 ; ii < nee ; ii++ ) free(ename[ii]) ;
207 
208      if( num_anlz > 0 ){                               /* got some! */
209        STATUS("opening ANALYZE files") ;
210        for( ii=0 ; ii < num_anlz ; ii++ ){             /* loop over files */
211          dset = THD_open_analyze( fn_anlz[ii] ) ;      /* try it on */
212          if( !ISVALID_DSET(dset) ) continue ;          /* doesn't fit? */
213          nds = sess->num_dsset ;
214          if( nds >= THD_MAX_SESSION_SIZE ){
215            fprintf(stderr,
216              "\n*** Session %s table overflow with ANALYZE dataset %s ***\n",
217              sessname , fn_anlz[ii] ) ;
218            THD_delete_3dim_dataset( dset , False ) ;
219            break ; /* out of for loop */
220         }
221         iview = dset->view_type ;
222         SET_SESSION_DSET(dset, sess, nds, iview);
223         sess->num_dsset ++ ;
224 #ifdef ALLOW_FSL_FEAT
225              if( strcmp(DSET_PREFIX(dset),"example_func.hdr") == 0 ) feat_exf = nds;
226         else if( strcmp(DSET_PREFIX(dset),"highres.hdr")      == 0 ) feat_hrs = nds;
227         else if( strcmp(DSET_PREFIX(dset),"standard.hdr")     == 0 ) feat_std = nds;
228 #endif
229        } /* end of loop over files */
230 
231        MCW_free_expand( num_anlz , fn_anlz ) ;  /* free file list */
232 
233        /* now read in linear mappings (.mat files) for FSL FEAT */
234 
235 #ifdef ALLOW_FSL_FEAT
236        if( feat_exf >= 0                        &&
237            sess->num_dsset - feat_nds_start > 0 &&
238            (feat_hrs >= 0 || feat_std >= 0)       ){  /* found FEAT files */
239 
240          THD_3dim_dataset *dset_exf=NULL, *dset_hrs=NULL, *dset_std=NULL ;
241          char fnam[THD_MAX_NAME] ;
242          FILE *fp ;
243          float a11,a12,a13,s1 , a21,a22,a23,s2 , a31,a32,a33,s3 ;
244          THD_warp *warp_exf_hrs=NULL , *warp_exf_std=NULL , *warp_std_hrs=NULL ;
245 
246                              dset_exf = GET_SESSION_DSET(sess, feat_exf, 0) ;  /* Must have this. */
247          if( feat_hrs >= 0 ) dset_hrs = GET_SESSION_DSET(sess, feat_hrs, 0) ;  /* And at least   */
248          if( feat_std >= 0 ) dset_std = GET_SESSION_DSET(sess, feat_std, 0) ;  /* one of these. */
249 
250 /*                              dset_exf = sess->dsset_xform_table[feat_exf][0] ;
251          if( feat_hrs >= 0 ) dset_hrs = sess->dsset_xform_table[feat_hrs][0] ;
252          if( feat_std >= 0 ) dset_std = sess->dsset_xform_table[feat_std][0] ; */
253 
254          /* try to read the warp from example_func (EPI) to standard */
255 
256          if( dset_std != NULL ){
257            strcpy(fnam,sess->sessname) ; strcat(fnam,"example_func2standard.mat") ;
258            fp = fopen(fnam,"r") ;
259            if( fp != NULL ){
260              ii = fscanf(fp,"%f%f%f%f %f%f%f%f %f%f%f%f" ,
261                          &a11,&a12,&a13,&s1 , &a21,&a22,&a23,&s2 ,
262                                               &a31,&a32,&a33,&s3   ) ;
263              if( ii == 12 )
264                warp_exf_std = AFNI_make_affwarp_12( a11,a12,a13,s1 ,
265                                                     a21,a22,a23,s2 ,
266                                                     a31,a32,a33,s3   ) ;
267              fclose(fp) ;
268            }
269 
270            /* 28 Aug 2002:
271                (i)  correct orientation of example_func
272                (ii) correct warp for non-DICOM order of coords in datasets      */
273 
274            if( warp_exf_std != NULL ){
275              THD_mat33 mmm,nnn ;
276              int   ix , iy , iz ;
277              float ax , ay , az ;
278 
279 #if 0
280 printf("warp_exf_std BEFORE:") ; DUMP_LMAP(warp_exf_std->rig_bod.warp) ;
281 #endif
282 
283 #undef FIX_EXF
284 #ifdef FIX_EXF
285              /* make matrix that transforms standard to example_func */
286 
287              LOAD_MAT(mmm,a11,a12,a13,a21,a22,a23,a31,a32,a33) ;
288              nnn = MAT_INV(mmm) ;
289              UNLOAD_MAT(nnn,a11,a12,a13,a21,a22,a23,a31,a32,a33) ;
290 
291              /* for each for, find index and value of largest element */
292 
293              ix = 1 ; ax = a11 ;
294              if( fabs(a12) > fabs(ax) ){ ix = 2 ; ax = a12 ; }
295              if( fabs(a13) > fabs(ax) ){ ix = 3 ; ax = a13 ; }
296 
297              iy = 1 ; ay = a21 ;
298              if( fabs(a22) > fabs(ay) ){ iy = 2 ; ay = a22 ; }
299              if( fabs(a23) > fabs(ay) ){ iy = 3 ; ay = a23 ; }
300 
301              iz = 1 ; az = a31 ;
302              if( fabs(a32) > fabs(az) ){ iz = 2 ; az = a32 ; }
303              if( fabs(a33) > fabs(az) ){ iz = 3 ; az = a33 ; }
304 #else
305              ix = 1 ; iy = 2 ; iz = 3 ;
306 #endif
307 
308              if( ix+iy+iz == 6 ){  /* (ix,iy,iz) must be a permutation of (1,2,3) */
309                THD_ivec3 orixyz ;
310                THD_fvec3 dxyz ;
311                THD_dataxes *daxes = dset_exf->daxes ;
312                THD_warp *from_exf , *to_std ;
313 
314 #ifdef FIX_EXF
315                /** fix orientation of dset_exf dataset **/
316 
317                switch(ix){
318                  case 1:    /* example_func x-axis is mainly same
319                                as standard x-axis (ax>0) or its opposite (ax<0) */
320                    if( ax > 0 ) ix = dset_std->daxes->xxorient ;
321                    else         ix = ORIENT_OPPOSITE(dset_std->daxes->xxorient) ;
322                  break ;
323 
324                  case 2:    /* example_func x-axis is mainly same
325                                as standard y-axis (ax>0) or its opposite (ax<0) */
326                    if( ax > 0 ) ix = dset_std->daxes->yyorient ;
327                    else         ix = ORIENT_OPPOSITE(dset_std->daxes->yyorient) ;
328                  break ;
329 
330                  case 3:    /* example_func x-axis is mainly same
331                                as standard z-axis (ax>0) or its opposite (ax<0) */
332                    if( ax > 0 ) ix = dset_std->daxes->zzorient ;
333                    else         ix = ORIENT_OPPOSITE(dset_std->daxes->zzorient) ;
334                  break ;
335                }
336                switch(iy){
337                  case 1: if( ay > 0 ) iy = dset_std->daxes->xxorient ;
338                          else         iy = ORIENT_OPPOSITE(dset_std->daxes->xxorient) ;
339                  break ;
340                  case 2: if( ay > 0 ) iy = dset_std->daxes->yyorient ;
341                          else         iy = ORIENT_OPPOSITE(dset_std->daxes->yyorient) ;
342                  break ;
343                  case 3: if( ay > 0 ) iy = dset_std->daxes->zzorient ;
344                          else         iy = ORIENT_OPPOSITE(dset_std->daxes->zzorient) ;
345                  break ;
346                }
347                switch(iz){
348                  case 1: if( az > 0 ) iz = dset_std->daxes->xxorient ;
349                          else         iz = ORIENT_OPPOSITE(dset_std->daxes->xxorient) ;
350                  break ;
351                  case 2: if( az > 0 ) iz = dset_std->daxes->yyorient ;
352                          else         iz = ORIENT_OPPOSITE(dset_std->daxes->yyorient) ;
353                  break ;
354                  case 3: if( az > 0 ) iz = dset_std->daxes->zzorient ;
355                          else         iz = ORIENT_OPPOSITE(dset_std->daxes->zzorient) ;
356                  break ;
357                }
358                orixyz.ijk[0] = ix ; orixyz.ijk[1] = iy ; orixyz.ijk[2] = iz ;
359                dxyz.xyz[0] = fabs(daxes->xxdel) ;
360                dxyz.xyz[1] = fabs(daxes->yydel) ;
361                dxyz.xyz[2] = fabs(daxes->zzdel) ;
362                if( ORIENT_sign[ix] == '-' ) dxyz.xyz[0] = -dxyz.xyz[0] ;
363                if( ORIENT_sign[iy] == '-' ) dxyz.xyz[1] = -dxyz.xyz[1] ;
364                if( ORIENT_sign[iz] == '-' ) dxyz.xyz[2] = -dxyz.xyz[2] ;
365                EDIT_dset_items( dset_exf , ADN_xyzorient,orixyz, ADN_xyzdel,dxyz, ADN_none ) ;
366 #endif
367 
368                /** now fix warp from example_func to standard to do it in DICOM coords **/
369 
370                nnn = SNGL_mat_to_dicomm( dset_exf ) ; mmm = TRANSPOSE_MAT(nnn) ;
371                from_exf = AFNI_make_affwarp_mat( mmm ) ;
372                nnn = SNGL_mat_to_dicomm( dset_std ) ;
373                to_std = AFNI_make_affwarp_mat( nnn ) ;
374                AFNI_concatenate_warp( warp_exf_std , from_exf ) ;
375                AFNI_concatenate_warp( to_std , warp_exf_std ) ;
376                myRwcFree(warp_exf_std) ; myRwcFree(from_exf) ;
377                warp_exf_std = to_std ;
378 
379 #if 0
380 printf("warp_exf_std AFTER:") ; DUMP_LMAP(warp_exf_std->rig_bod.warp) ;
381 #endif
382 
383              } /* end of if warp had reasonable axis combinations */
384            } /* end of if we got a good warp */
385          }
386 
387          /* try to read the warp from example_func (EPI) to highres */
388 
389          if( dset_hrs != NULL ){
390            strcpy(fnam,sess->sessname) ; strcat(fnam,"example_func2highres.mat") ;
391            fp = fopen(fnam,"r") ;
392            if( fp != NULL ){
393              ii = fscanf(fp,"%f%f%f%f %f%f%f%f %f%f%f%f" ,
394                          &a11,&a12,&a13,&s1 , &a21,&a22,&a23,&s2 ,
395                                               &a31,&a32,&a33,&s3   ) ;
396              if( ii == 12 )
397                warp_exf_hrs = AFNI_make_affwarp_12( a11,a12,a13,s1 ,
398                                                     a21,a22,a23,s2 ,
399                                                     a31,a32,a33,s3   ) ;
400              fclose(fp) ;
401            }
402 
403            /* 28 Aug 2002: correct warp for non-DICOM order of coords in datasets */
404 
405            if( warp_exf_hrs != NULL ){
406              THD_mat33 mmm,nnn ;
407              THD_warp *from_exf , *to_hrs ;
408 
409 #if 0
410 printf("warp_exf_hrs BEFORE:") ; DUMP_LMAP(warp_exf_hrs->rig_bod.warp) ;
411 #endif
412 
413              nnn = SNGL_mat_to_dicomm( dset_exf ) ; mmm = TRANSPOSE_MAT(nnn) ;
414              from_exf = AFNI_make_affwarp_mat( mmm ) ;
415              nnn = SNGL_mat_to_dicomm( dset_hrs ) ;
416              to_hrs = AFNI_make_affwarp_mat( nnn ) ;
417              AFNI_concatenate_warp( warp_exf_hrs , from_exf ) ;
418              AFNI_concatenate_warp( to_hrs , warp_exf_hrs ) ;
419              myRwcFree(warp_exf_hrs) ; myRwcFree(from_exf) ;
420              warp_exf_hrs = to_hrs ;
421 
422 #if 0
423 printf("warp_exf_hrs AFTER:") ; DUMP_LMAP(warp_exf_hrs->rig_bod.warp) ;
424 #endif
425            }
426          }
427 
428          /* try to read the warp from standard to highres */
429 
430          if( dset_hrs != NULL && dset_std != NULL ){
431            strcpy(fnam,sess->sessname) ; strcat(fnam,"standard2highres.mat") ;
432            fp = fopen(fnam,"r") ;
433            if( fp != NULL ){
434              ii = fscanf(fp,"%f%f%f%f %f%f%f%f %f%f%f%f" ,
435                          &a11,&a12,&a13,&s1 , &a21,&a22,&a23,&s2 ,
436                                               &a31,&a32,&a33,&s3   ) ;
437              if( ii == 12 )
438                warp_std_hrs = AFNI_make_affwarp_12( a11,a12,a13,s1 ,
439                                                     a21,a22,a23,s2 ,
440                                                     a31,a32,a33,s3   ) ;
441              fclose(fp) ;
442            }
443 
444            /* 28 Aug 2002: correct warp for non-DICOM order of coords in datasets */
445 
446            if( warp_std_hrs != NULL ){
447              THD_mat33 mmm,nnn ;
448              THD_warp *from_std , *to_hrs ;
449 
450 #if 0
451 printf("warp_std_hrs BEFORE:") ; DUMP_LMAP(warp_std_hrs->rig_bod.warp) ;
452 #endif
453 
454              nnn = SNGL_mat_to_dicomm( dset_std ) ; mmm = TRANSPOSE_MAT(nnn) ;
455              from_std = AFNI_make_affwarp_mat( mmm ) ;
456              nnn = SNGL_mat_to_dicomm( dset_hrs ) ;
457              to_hrs = AFNI_make_affwarp_mat( nnn ) ;
458              AFNI_concatenate_warp( warp_std_hrs , from_std ) ;
459              AFNI_concatenate_warp( to_hrs , warp_std_hrs ) ;
460              myRwcFree(warp_std_hrs) ; myRwcFree(from_std) ;
461              warp_std_hrs = to_hrs ;
462 
463 #if 0
464 printf("warp_std_hrs AFTER:") ; DUMP_LMAP(warp_std_hrs->rig_bod.warp) ;
465 #endif
466            }
467          }
468 
469          /* if we have these warps,
470             then build a hashtable describing their use in
471             transforming funcs (in exf coords) to hrs and std coords */
472 
473          if( warp_exf_hrs != NULL || warp_exf_std != NULL ){
474 
475            if( sess->warptable == NULL )
476               sess->warptable = new_Htable(0) ;  /* use minimum table size */
477 
478            if( warp_exf_hrs != NULL ){
479              for( ii=feat_nds_start ; ii < sess->num_dsset ; ii++ ){
480                temp_dset = GET_SESSION_DSET(sess, ii, 0);
481                if( ISFUNC(temp_dset) ){
482 /*               if( ISFUNC(sess->dsset_xform_table[ii][0]) ){*/
483 /*                 sprintf(fnam,"%s,%s",dset_hrs->idcode.str,sess->dsset_xform_table[ii][0]->idcode.str) ;*/
484                  sprintf(fnam,"%s,%s",dset_hrs->idcode.str,temp_dset->idcode.str) ;
485                  addto_Htable( fnam , warp_exf_hrs , sess->warptable ) ;
486                }
487              }
488              for( ii=feat_nds_start ; ii < sess->num_dsset ; ii++ ){
489                temp_dset = GET_SESSION_DSET(sess, ii, 0);
490                if( ii != feat_hrs && ii != feat_std && ISANAT(temp_dset) ){
491 /*               if( ii != feat_hrs && ii != feat_std && ISANAT(sess->dsset_xform_table[ii][0]) ){*/
492 /*                 sprintf(fnam,"%s,%s",dset_hrs->idcode.str,sess->dsset_xform_table[ii][0]->idcode.str) ;*/
493                  sprintf(fnam,"%s,%s",dset_hrs->idcode.str,temp_dset->idcode.str) ;
494                  addto_Htable( fnam , warp_exf_hrs , sess->warptable ) ;
495                }
496              }
497            }
498 
499            if( warp_exf_std != NULL ){
500              for( ii=feat_nds_start ; ii < sess->num_dsset ; ii++ ){
501                temp_dset = GET_SESSION_DSET(sess, ii, 0);
502                if( ISFUNC(temp_dset) ){
503 /*             if( ISFUNC(sess->dsset_xform_table[ii][0]) ){*/
504 /*                sprintf(fnam,"%s,%s",dset_std->idcode.str,sess->dsset_xform_table[ii][0]->idcode.str) ;*/
505                  sprintf(fnam,"%s,%s",dset_std->idcode.str,temp_dset->idcode.str) ;
506                  addto_Htable( fnam , warp_exf_std , sess->warptable ) ;
507                }
508              }
509              for( ii=feat_nds_start ; ii < sess->num_dsset ; ii++ ){
510                temp_dset = GET_SESSION_DSET(sess, ii, 0);
511                if( ii != feat_hrs && ii != feat_std && ISANAT(temp_dset) ){
512                  sprintf(fnam,"%s,%s",dset_std->idcode.str,temp_dset->idcode.str) ;
513 /*               if( ii != feat_hrs && ii != feat_std && ISANAT(sess->dsset_xform_table[ii][0]) ){ */
514 /*                 sprintf(fnam,"%s,%s",dset_std->idcode.str,sess->dsset_xform_table[ii][0]->idcode.str) ; */
515                  addto_Htable( fnam , warp_exf_std , sess->warptable ) ;
516                }
517              }
518            }
519 
520            if( warp_std_hrs != NULL ){
521              sprintf(fnam,"%s,%s",dset_hrs->idcode.str,dset_std->idcode.str) ;
522              addto_Htable( fnam , warp_std_hrs , sess->warptable ) ;
523            }
524 
525          } /* end of making warptable */
526        } /* end of FEATing */
527 #endif
528 
529      } /* end of if we found ANALYZE files */
530    }
531 
532    /*-- 04 Dec 2002: try to read CTF .mri and .svl "datasets" --*/
533 
534    if( !AFNI_noenv("AFNI_CTF_DATASETS") ){
535      char *ename[2] , **fn_ctf ;
536      int num_ctf , ii ;
537 
538      STATUS("looking for CTF files") ;
539 
540      ename[0] = AFMALL(char, THD_MAX_NAME) ;
541      ename[1] = AFMALL(char, THD_MAX_NAME) ;
542      strcpy(ename[0],sess->sessname) ; strcat(ename[0],"*.mri") ;
543      strcpy(ename[1],sess->sessname) ; strcat(ename[1],"*.svl") ;
544      MCW_file_expand( 2,ename , &num_ctf,&fn_ctf ) ;  /* find files */
545      free(ename[0]) ; free(ename[1]) ;
546 
547      if( num_ctf > 0 ){                               /* got some files! */
548        STATUS("opening CTF files") ;
549        for( ii=0 ; ii < num_ctf ; ii++ ){             /* loop over files */
550 
551          if( strstr(fn_ctf[ii],".mri") != NULL )      /* try to read: */
552            dset = THD_open_ctfmri( fn_ctf[ii] ) ;     /*   as MRI */
553          else if( strstr(fn_ctf[ii],".svl") != NULL )
554            dset = THD_open_ctfsam( fn_ctf[ii] ) ;     /*   as SAM */
555 
556          if( !ISVALID_DSET(dset) ) continue ;         /* doesn't read? */
557          nds = sess->num_dsset ;
558          if( nds >= THD_MAX_SESSION_SIZE ){
559            fprintf(stderr,
560             "\n*** Session %s table overflow with dataset %s ***\n",
561             sessname , fn_ctf[ii] ) ;
562            THD_delete_3dim_dataset( dset , False ) ;
563            break ; /* out of for loop */
564          }
565          iview = dset->view_type ;
566          SET_SESSION_DSET(dset, sess, nds, iview);
567          sess->num_dsset ++ ;
568        } /* end of loop over files */
569        MCW_free_expand( num_ctf , fn_ctf ) ;
570      } /* end of if we found CTF files */
571    }
572 
573    /*-- 02 Feb 2018: try to read .jpg and .png "datasets" --*/
574 
575    if( !AFNI_noenv("AFNI_IMAGE_DATASETS") ){
576      char *ename[4] , **fn_img ;
577      int num_img , ii ;
578 
579      STATUS("looking for JPG/PNG files") ;
580 
581      ename[0] = AFMALL(char, THD_MAX_NAME) ;
582      ename[1] = AFMALL(char, THD_MAX_NAME) ;
583      ename[2] = AFMALL(char, THD_MAX_NAME) ;
584      ename[3] = AFMALL(char, THD_MAX_NAME) ;
585      strcpy(ename[0],sess->sessname) ; strcat(ename[0],"*.png") ;
586      strcpy(ename[1],sess->sessname) ; strcat(ename[1],"*.PNG") ;
587      strcpy(ename[2],sess->sessname) ; strcat(ename[2],"*.jpg") ;
588      strcpy(ename[3],sess->sessname) ; strcat(ename[3],"*.JPG") ;
589      MCW_file_expand( 4,ename , &num_img,&fn_img ) ;  /* find files */
590      free(ename[0]) ; free(ename[1]) ; free(ename[2]) ; free(ename[3]) ;
591 
592      if( num_img > 0 ){                               /* got some files! */
593        STATUS("opening JPG/PNG files") ;
594        for( ii=0 ; ii < num_img ; ii++ ){             /* loop over files */
595 
596          dset = THD_open_image( fn_img[ii] ) ;        /* try top read */
597 
598          if( !ISVALID_DSET(dset) ) continue ;         /* doesn't read? */
599          nds = sess->num_dsset ;
600          if( nds >= THD_MAX_SESSION_SIZE ){
601            fprintf(stderr,
602             "\n*** Session %s table overflow with dataset %s ***\n",
603             sessname , fn_img[ii] ) ;
604            THD_delete_3dim_dataset( dset , False ) ;
605            break ; /* out of for(ii) loop */
606          }
607          iview = dset->view_type ;
608          SET_SESSION_DSET(dset, sess, nds, iview);
609          sess->num_dsset ++ ;
610        } /* end of loop over files */
611        MCW_free_expand( num_img , fn_img ) ;
612      } /* end of if we found JPG/PNG files */
613    }
614 
615    /*-- 03 Dec 2001: try to read MPEG "datasets" --*/
616 
617    if( AFNI_yesenv("AFNI_MPEG_DATASETS") ){
618      char ename[4*THD_MAX_NAME+64] , **fn_mpeg ;
619      int num_mpeg , ii ;
620 
621      STATUS("looking for MPEG files") ;
622 
623      sprintf(ename,"%s*.mpg %s*.mpeg %s*.MPEG %s*.MPG" ,
624              sess->sessname, sess->sessname, sess->sessname, sess->sessname ) ;
625      MCW_wildcards( ename , &num_mpeg , &fn_mpeg ) ;   /* find files */
626 
627      if( num_mpeg > 0 ){                               /* got some! */
628        STATUS("opening MPEG files") ;
629        for( ii=0 ; ii < num_mpeg ; ii++ ){             /* loop over files */
630          dset = THD_open_mpeg( fn_mpeg[ii] ) ;         /* try it on */
631          if( !ISVALID_DSET(dset) ) continue ;          /* doesn't fit? */
632          nds = sess->num_dsset ;
633          if( nds >= THD_MAX_SESSION_SIZE ){
634            fprintf(stderr,
635              "\n*** Session %s table overflow with MPEG dataset %s ***\n",
636              sessname , fn_mpeg[ii] ) ;
637            THD_delete_3dim_dataset( dset , False ) ;
638            break ; /* out of for loop */
639          }
640          iview = dset->view_type ;
641          SET_SESSION_DSET(dset, sess, nds, iview);
642          sess->num_dsset ++ ;
643        } /* end of loop over files */
644        MCW_free_expand( num_mpeg , fn_mpeg ) ;
645      } /* end of if we found MPEG files */
646    }
647 
648    /*-- done! --*/
649 
650    if( sess->num_dsset == 0 ){
651       myRwcFree( sess ) ; RETURN( NULL ) ;
652    }
653 
654    /*-- 29 Jul 2003: order dataset rows --*/
655 
656    THD_order_session( sess ) ;
657 
658    RETURN( sess ) ;
659 }
660 
661 /*-------------------------------------------------------------------------*/
662 /*! Order the datasets within a session (anats first, funcs last).
663 ---------------------------------------------------------------------------*/
664 
THD_order_session(THD_session * sess)665 void THD_order_session( THD_session *sess )
666 {
667    THD_3dim_dataset *qset[THD_MAX_SESSION_SIZE][LAST_VIEW_TYPE+1] ;
668    THD_3dim_dataset *dset=NULL ;
669    int iview , ids , nds ;
670 
671 ENTRY("THD_order_session") ;
672    if( sess == NULL || sess->num_dsset <= 1 ) EXRETURN ;
673 
674    /* put anats into qset */
675 
676    nds = 0 ;
677    for( ids=0 ; ids < sess->num_dsset ; ids++ ){
678      for( iview=0 ; iview <= LAST_VIEW_TYPE ; iview++ ){
679        dset = GET_SESSION_DSET(sess, ids, iview);
680 /*       dset = sess->dsset_xform_table[ids][iview] ;*/
681        if( dset != NULL && ISANAT(dset) ) break ;
682      }
683      if( iview <= LAST_VIEW_TYPE ){
684        for( iview=0 ; iview <= LAST_VIEW_TYPE ; iview++ )
685          qset[nds][iview] = GET_SESSION_DSET(sess, ids, iview);
686 /*       qset[nds][iview] = sess->dsset_xform_table[ids][iview] ;*/
687        nds++ ;
688      }
689    }
690 
691    /* put funcs into qset */
692 
693    for( ids=0 ; ids < sess->num_dsset ; ids++ ){
694      for( iview=0 ; iview <= LAST_VIEW_TYPE ; iview++ ){
695        dset = GET_SESSION_DSET(sess, ids, iview);
696 /*     dset = sess->dsset_xform_table[ids][iview] ;*/
697        /* if( dset != NULL && ISFUNC(dset) ) break ; */
698        if( dset ) break ;
699      }
700 
701      /* might be ANAT in one view and FUNC in another (consider qsort()?)
702       *                                        17 Nov 2011 [rickr, dglen] */
703      if( ! dset )            continue;
704      else if( ISANAT(dset) ) continue;   /* already found above */
705      /* else, add ... */
706 
707      if( iview <= LAST_VIEW_TYPE ){
708        for( iview=0 ; iview <= LAST_VIEW_TYPE ; iview++ )
709           qset[nds][iview] = GET_SESSION_DSET(sess, ids, iview);
710 /*       qset[nds][iview] = sess->dsset_xform_table[ids][iview] ;*/
711        nds++ ;
712      }
713    }
714 
715    /* copy qset back into sess, overwriting dsset */
716 
717    for( ids=0 ; ids < nds ; ids++ )
718      for( iview=0 ; iview <= LAST_VIEW_TYPE ; iview++ )
719         SET_SESSION_DSET(qset[ids][iview], sess, ids, iview);
720    sess->num_dsset = nds ;  /* shouldn't change */
721    EXRETURN ;
722 }
723 
724 /*---------------------------------------------------------------------*/
725 /*! Append datasets in THD_session ssb to those in ssa.
726     \date 20 Dec 2001
727 */
728 
THD_append_sessions(THD_session * ssa,THD_session * ssb)729 void THD_append_sessions( THD_session *ssa , THD_session *ssb )
730 {
731    int qs, qd, vv ;
732    THD_3dim_dataset *temp_dset;
733 
734 ENTRY("THD_append_sessions") ;
735 
736    if( !ISVALID_SESSION(ssa) || !ISVALID_SESSION(ssb)  ) EXRETURN ;
737    if( THD_equiv_files(ssa->sessname,ssb->sessname)==1 ) EXRETURN ;
738 
739    qs = ssa->num_dsset ;
740    for( qd=0; qd < ssb->num_dsset && qd+qs < THD_MAX_SESSION_SIZE ; qd++ )
741      for( vv=0 ; vv <= LAST_VIEW_TYPE ; vv++ ) {
742        temp_dset = GET_SESSION_DSET(ssb, qd, vv);
743        SET_SESSION_DSET(temp_dset, ssa,qd+qs,vv);
744      }
745    ssa->num_dsset += qd ;
746 
747    EXRETURN ;
748 }
749 
750 /*---------------------------------------------------------------------*/
751 /* read from the pipe into 1 big string [01 Feb 2018] */
752 
753 #undef  PSIZ
754 #define PSIZ 8192
755 
THD_suck_pipe(char * cmd)756 char * THD_suck_pipe( char *cmd )
757 {
758    char *flist=NULL ; FILE *fp=NULL ;
759 
760    if( cmd == NULL || *cmd == '\0' ) return NULL ;
761 
762    fp = popen( cmd , "r" ) ;
763    if( fp == NULL ){
764      ERROR_message("Unable to pipe command '%s'",cmd) ;
765      return NULL ;
766    }
767 
768    flist = (char *)malloc(sizeof(char)*PSIZ) ; flist[0] = '\0' ;
769    while( fgets(flist+strlen(flist),PSIZ-2,fp) != NULL ){
770      flist = (char *)realloc(flist,sizeof(char)*(strlen(flist)+PSIZ)) ;
771    }
772 
773    (void)pclose(fp) ;
774    if( *flist == '\0' ){ free(flist); flist=NULL; }
775 
776    return flist ;
777 }
778 
779 /*---------------------------------------------------------------------*/
780 /* Find all sub-directories with BIDS name form sub-XXX [01 Feb 2018] */
781 
THD_get_subdirs_bysub(char * dirname,char * subid)782 NI_str_array * THD_get_subdirs_bysub( char *dirname , char *subid )
783 {
784    NI_str_array *qsar=NULL ;
785    char *cmd , *flist , *newsubid=NULL ; FILE *fp ;
786 
787 ENTRY("THD_get_subdirs_bysub") ;
788 
789    if( ! THD_is_directory(dirname) || subid == NULL || strlen(subid) < 2 )
790      RETURN(NULL) ;
791 
792    if( strncmp(subid,"sub-",4) != 0 ){
793      newsubid = (char *)malloc(sizeof(char)*(8+strlen(subid))) ;
794      sprintf(newsubid,"sub-%s",subid) ;
795    } else {
796      newsubid = strdup(subid) ;
797    }
798 
799    /* get dirnames via piping from find */
800 
801    cmd = (char *)malloc(sizeof(char)*(128+strlen(dirname)+strlen(newsubid))) ;
802    sprintf( cmd, "find %s -maxdepth 9 -type d -name '%s'", dirname , newsubid ) ;
803 
804    flist = THD_suck_pipe( cmd ) ;
805 
806    if( flist == NULL || strlen(flist) < 4 ){
807 #if 0
808      ERROR_message("-bysub: Didn't find any '%s' under %s",newsubid,dirname) ;
809 #endif
810      free(newsubid) ; free(cmd) ;
811      if( flist != NULL ) free(flist) ;
812      RETURN(NULL) ;
813    }
814 
815    /* break output into discrete strings */
816 
817    qsar = NI_decode_string_list( flist , ";" ) ;
818    if( qsar == NULL || qsar->num == 0 ){    /* should never happen */
819      ERROR_message("-bysub: Unable to decode output from command \"%s\"",cmd) ;
820      free(newsubid) ; free(cmd) ; free(flist) ; RETURN(NULL) ;
821    }
822 
823 #if 0
824   { int qq ;
825     INFO_message("-bysub: Found %d '%s' under %s",qsar->num,newsubid,dirname) ;
826     for( qq=0 ; qq < qsar->num ; qq++ ) ININFO_message(" %s",qsar->str[qq]) ;
827   }
828 #endif
829 
830    free(newsubid) ; free(cmd) ; free(flist) ; RETURN(qsar) ;
831 }
832 
833 /*---------------------------------------------------------------------*/
834 /* A recursive form of THD_init_session() [01 Feb 2018] */
835 
THD_init_session_recursive(char * dirname)836 THD_session * THD_init_session_recursive( char *dirname )
837 {
838    NI_str_array *qsar ;
839    THD_session *sess=NULL , *qss ;
840    char *cmd , *flist ;
841    int qq , nss=0 ;
842 ENTRY("THD_init_session_recursive") ;
843 
844    if( dirname == NULL || *dirname == '\0' ) RETURN(NULL) ;
845 
846    /* get list of all subdir names */
847 
848    cmd = (char *)malloc(sizeof(char)*(strlen(dirname)+128)) ;
849    sprintf( cmd, "find %s -type d -depth -9", dirname) ;
850    flist = THD_suck_pipe( cmd ) ;
851    if( flist == NULL || strlen(flist) < 1 ) RETURN(NULL) ;
852 
853    /* break output into discrete strings */
854 
855    qsar = NI_decode_string_list( flist , ";" ) ;
856    free(flist) ;
857    if( qsar == NULL || qsar->num == 0 ){    /* should never happen */
858      ERROR_message("Unable to decode output from command \"%s\"",cmd) ;
859      free(cmd) ; RETURN(NULL) ;
860    }
861 
862    /* try to read each subdir as a session by itself */
863 
864    for( qq=0 ; qq < qsar->num ; qq++ ){
865      if( qsar->str[qq] != NULL && qsar->str[qq][0] != '\0' ){
866        qss = THD_init_session( qsar->str[qq] ) ;
867        if( qss != NULL ){    /* got something */
868          if( sess == NULL ){ /* the first something */
869            sess = qss ;
870          } else {            /* a later something */
871            THD_append_sessions( sess , qss ) ;
872          }
873          nss++ ;
874        }
875      }
876    }
877 
878    if( nss > 1 ){
879      THD_order_session( sess ) ;
880      sess->is_collection = 1 ;
881    }
882 
883    free(cmd) ;
884    RETURN(sess) ;
885 }
886 
887 /*---------------------------------------------------------------------*/
888 /* Make a session from all subdirectories that match the subject ID */
889 
THD_init_session_bysub(char * dirname,char * subid)890 THD_session * THD_init_session_bysub( char *dirname , char *subid )
891 {
892    NI_str_array *qsar ;
893    THD_session *sess=NULL , *qss ;
894    int qq , nss=0 ;
895 
896 ENTRY("THD_init_session_bysub") ;
897 
898    qsar = THD_get_subdirs_bysub( dirname , subid ) ;
899    if( qsar == NULL ) RETURN(NULL) ;
900 
901    for( qq=0 ; qq < qsar->num ; qq++ ){
902      if( qsar->str[qq] != NULL && qsar->str[qq][0] != '\0' ){
903        qss = THD_init_session_recursive( qsar->str[qq] ) ;
904        if( qss != NULL ){    /* got something */
905          if( sess == NULL ){ /* the first something */
906            sess = qss ;
907          } else {            /* a later something */
908            THD_append_sessions( sess , qss ) ;
909          }
910          nss++ ;
911        }
912      }
913    }
914    NI_delete_str_array(qsar) ;  /* oops [15 Apr 2019] */
915 
916    if( sess != NULL && nss > 1 ){
917      THD_order_session( sess ) ;
918      sess->is_collection = 1 ;
919    }
920    if( sess != NULL ){
921      char *env = my_getenv( "AFNI_SESSTRAIL" ) ; int tt = 0 ;
922      if( env != NULL ) tt = (int)strtol(env,NULL,10) ;
923      env = THD_trailname(dirname,tt) ;
924      sprintf( sess->lastname , "%.32s:%s",env,subid) ;
925    }
926 
927    RETURN(sess) ;
928 }
929