00001 #include "mrilib.h"
00002 #include "afni.h"
00003 
00004 static int           have_dseTT = -1   ;
00005 static THD_3dim_dataset * dseTT = NULL ;
00006 static THD_3dim_dataset * dseTT_big = NULL ; /* 01 Aug 2001 */
00007 
00008 #define MAX_FIND 9                    /* max number to find within WAMIRAD  */
00009 #define WAMIRAD  7.5                  /* search radius: must not exceed 9.5 */
00010 static MCW_cluster * wamiclust=NULL ;
00011 
00012 /*-----------------------------------------------------------------------*/
00013 
00014 THD_3dim_dataset * TT_retrieve_atlas(void)
00015 {
00016    if( have_dseTT < 0 ) TT_load_atlas() ;
00017    return dseTT ;                         /* might be NULL */
00018 }
00019 
00020 /*-----------------------------------------------------------------------*/
00021 
00022 THD_3dim_dataset * TT_retrieve_atlas_big(void) /* 01 Aug 2001 */
00023 {
00024    if( dseTT_big != NULL ) return dseTT_big ;
00025    if( have_dseTT < 0    ) TT_load_atlas() ;
00026    if( dseTT == NULL     ) return NULL ;
00027    dseTT_big = THD_zeropad( dseTT , 10,0,0,0,0,0 , "TTatlas_big" , 0 ) ;
00028    DSET_unload( dseTT ) ; /* probably won't need again */
00029    return dseTT_big ;
00030 }
00031 
00032 /*-----------------------------------------------------------------------*/
00033 
00034 THD_3dim_dataset * TT_retrieve_atlas_either(void) /* 22 Aug 2001 */
00035 {
00036    if( dseTT_big != NULL ) return dseTT_big ;
00037    if( dseTT     != NULL ) return dseTT     ;
00038    if( have_dseTT < 0    ) TT_load_atlas()  ;
00039    return dseTT ;
00040 }
00041 
00042 /*-----------------------------------------------------------------------*/
00043 
00044 int TT_load_atlas(void)
00045 {
00046    char *epath, *elocal, ename[THD_MAX_NAME], dname[THD_MAX_NAME], *eee ;
00047    int epos , ll , ii , id ;
00048 
00049 ENTRY("TT_load_atlas") ;
00050 
00051    if( have_dseTT >= 0 ) RETURN(have_dseTT) ;  /* for later calls */
00052 
00053    have_dseTT = 0 ;  /* don't have it yet */
00054 
00055    /*----- 20 Aug 2001: see if user specified alternate database -----*/
00056 
00057    epath = getenv("AFNI_TTATLAS_DATASET") ;
00058    if( epath != NULL ){
00059      dseTT = THD_open_one_dataset( epath ) ;  /* try to open it */
00060      if( dseTT != NULL ){                     /* got it!!! */
00061        have_dseTT = 1; RETURN(1);
00062      }
00063    }
00064 
00065    /*----- get path to search -----*/
00066 
00067                        epath = getenv("AFNI_PLUGINPATH") ;
00068    if( epath == NULL ) epath = getenv("AFNI_PLUGIN_PATH") ;
00069    if( epath == NULL ) epath = getenv("PATH") ;
00070    if( epath == NULL ) RETURN(0) ;
00071 
00072    /*----- copy path list into local memory -----*/
00073 
00074    ll = strlen(epath) ;
00075    elocal = AFMALL(char, sizeof(char) * (ll+2) ) ;
00076 
00077    /*----- put a blank at the end -----*/
00078 
00079    strcpy( elocal , epath ) ; elocal[ll] = ' ' ; elocal[ll+1] = '\0' ;
00080 
00081    /*----- replace colons with blanks -----*/
00082 
00083    for( ii=0 ; ii < ll ; ii++ )
00084      if( elocal[ii] == ':' ) elocal[ii] = ' ' ;
00085 
00086    /*----- extract blank delimited strings;
00087            use as directory names to look for atlas -----*/
00088 
00089    epos = 0 ;
00090 
00091    do{
00092       ii = sscanf( elocal+epos , "%s%n" , ename , &id ); /* next substring */
00093       if( ii < 1 ) break ;                               /* none -> done   */
00094 
00095       epos += id ;                                 /* char after last scanned */
00096 
00097       ii = strlen(ename) ;                         /* make sure name has   */
00098       if( ename[ii-1] != '/' ){                    /* a trailing '/' on it */
00099           ename[ii]  = '/' ; ename[ii+1] = '\0' ;
00100       }
00101       strcpy(dname,ename) ;
00102       strcat(dname,"TTatlas+tlrc") ;               /* add dataset name */
00103 
00104       dseTT = THD_open_one_dataset( dname ) ;      /* try to open it */
00105 
00106       if( dseTT != NULL ){                         /* got it!!! */
00107          free(elocal); have_dseTT = 1; RETURN(1);
00108       }
00109 
00110       strcpy(dname,ename) ;           /* 21 Jul 2005: NIfTI-1.1 version? */
00111       strcat(dname,"TTatlas.nii.gz") ;
00112       dseTT = THD_open_one_dataset( dname ) ;
00113       if( dseTT != NULL ){ free(elocal); have_dseTT = 1; RETURN(1); }
00114 
00115    } while( epos < ll ) ;  /* scan until 'epos' is after end of epath */
00116 
00117    free(elocal) ; RETURN(0) ; /* got here -> didn't find it */
00118 }
00119 
00120 /*----------------------------------------------------------------------
00121   Allows the program to purge the memory used by the TT atlas dataset
00122 ------------------------------------------------------------------------*/
00123 
00124 void TT_purge_atlas(void)
00125 {
00126   PURGE_DSET(dseTT) ; return ;
00127 }
00128 
00129 void TT_purge_atlas_big(void)
00130 {
00131    if( dseTT_big != NULL ){ DSET_delete(dseTT_big) ; dseTT_big = NULL ; }
00132    return ;
00133 }
00134 
00135 /*----------------------------------------------------------------------
00136    Return a multi-line string of TT atlas labels near the given point
00137    (xx,yy,zz are in Dicom order coordinates).
00138    If NULL is returned, an error happened.  If no labels are near the
00139    given point, then a single line saying "you're lost" is returned.
00140    The string returned is malloc()-ed and should be free()-ed later.
00141    The string will end with a newline '\n' character.
00142 ------------------------------------------------------------------------*/
00143 
00144 char * TT_whereami( float xx , float yy , float zz )
00145 {
00146    int ii,kk , ix,jy,kz , nx,ny,nz,nxy , aa,bb,cc , ff,b2f,b4f,rff ;
00147    THD_ivec3 ijk ;
00148    byte *b2 , *b4 ;
00149    THD_string_array *sar ;
00150    char *b2lab , *b4lab ;
00151    char lbuf[256] , *rbuf ;
00152    int nfind, b2_find[MAX_FIND], b4_find[MAX_FIND], rr_find[MAX_FIND] ;
00153 
00154    THD_3dim_dataset * dset ; /* 01 Aug 2001 */
00155 
00156 ENTRY("TT_whereami") ;
00157 
00158    /*-- setup stuff: load atlas dataset, prepare search mask --*/
00159 
00160    if( dseTT == NULL ){
00161       ii = TT_load_atlas() ; if( ii == 0 ) RETURN(NULL) ;
00162    }
00163 
00164    /* 01 Aug 2001: maybe use big dataset (so don't need both in memory) */
00165 
00166    dset = (dseTT_big != NULL) ? dseTT_big : dseTT ;
00167 
00168 #if 0
00169 if( dset == dseTT_big ) fprintf(stderr,"TT_whereami using dseTT_big\n") ;
00170 else                    fprintf(stderr,"TT_whereami using dseTT\n") ;
00171 #endif
00172 
00173    DSET_load(dset) ;
00174    b2 = DSET_BRICK_ARRAY(dset,0) ; if( b2 == NULL ) RETURN(NULL) ;
00175    b4 = DSET_BRICK_ARRAY(dset,1) ; if( b4 == NULL ) RETURN(NULL) ;
00176 
00177    if( wamiclust == NULL ){
00178       wamiclust = MCW_build_mask( 0,0,0 , 1.0,1.0,1.0 , WAMIRAD ) ;
00179       if( wamiclust == NULL ) RETURN(NULL) ;  /* should not happen! */
00180 
00181       for( ii=0 ; ii < wamiclust->num_pt ; ii++ )       /* load radius */
00182          wamiclust->mag[ii] = (int)rint(sqrt((double)(
00183                                          wamiclust->i[ii]*wamiclust->i[ii]
00184                                         +wamiclust->j[ii]*wamiclust->j[ii]
00185                                         +wamiclust->k[ii]*wamiclust->k[ii]))) ;
00186 
00187       MCW_sort_cluster( wamiclust ) ;  /* sort by radius */
00188    }
00189 
00190    /*-- find locations near the given one that are in the Atlas --*/
00191 
00192    ijk = THD_3dmm_to_3dind( dset , TEMP_FVEC3(xx,yy,zz) ) ;  /* get indexes */
00193    UNLOAD_IVEC3(ijk,ix,jy,kz) ;                               /* from coords */
00194 
00195    nx = DSET_NX(dset) ;               /* size of TT atlas dataset axes */
00196    ny = DSET_NY(dset) ;
00197    nz = DSET_NZ(dset) ; nxy = nx*ny ;
00198 
00199    nfind = 0 ;
00200 
00201    /*-- check the exact input location --*/
00202 
00203    kk = ix + jy*nx + kz*nxy ;        /* index into brick arrays */
00204    if( b2[kk] != 0 || b4[kk] != 0 ){
00205       b2_find[0] = b2[kk] ;
00206       b4_find[0] = b4[kk] ;
00207       rr_find[0] = 0      ; nfind++ ;
00208    }
00209 
00210    /*-- check locations near it --*/
00211 
00212    for( ii=0 ; ii < wamiclust->num_pt ; ii++ ){
00213 
00214       /* compute index of nearby location, skipping if outside atlas */
00215 
00216       aa = ix + wamiclust->i[ii] ; if( aa < 0 || aa >= nx ) continue ;
00217       bb = jy + wamiclust->j[ii] ; if( bb < 0 || bb >= ny ) continue ;
00218       cc = kz + wamiclust->k[ii] ; if( cc < 0 || cc >= nz ) continue ;
00219 
00220       kk  = aa + bb*nx + cc*nxy ;   /* index into bricks */
00221       b2f = b2[kk] ; b4f = b4[kk] ; /* TT structures markers there */
00222 
00223       if( b2f == 0 && b4f == 0 )                            continue ;
00224 
00225       for( ff=0 ; ff < nfind ; ff++ ){       /* cast out         */
00226          if( b2f == b2_find[ff] ) b2f = 0 ;  /* duplicate labels */
00227          if( b4f == b4_find[ff] ) b4f = 0 ;  /* we already found */
00228       }
00229       if( b2f == 0 && b4f == 0 )                            continue ;
00230 
00231       b2_find[nfind] = b2f ;  /* save what we found */
00232       b4_find[nfind] = b4f ;
00233       rr_find[nfind] = (int) wamiclust->mag[ii] ;
00234       nfind++ ;
00235 
00236       if( nfind == MAX_FIND ) break ;  /* don't find TOO much */
00237    }
00238 
00239    /*-- assemble output string(s) --*/
00240 
00241 #define WAMI_HEAD "+++++++ nearby Talairach Daemon structures +++++++\n"
00242 #define WAMI_TAIL "\n******** Please use results with caution! ********"  \
00243                   "\n******** Brain anatomy is quite variable! ********"  \
00244                   "\n******** The database may contain errors! ********"
00245 
00246    if( nfind == 0 ){
00247       char xlab[24], ylab[24] , zlab[24] ;
00248       THD_fvec3 tv , mv ;
00249       float mx,my,mz ;
00250       char mxlab[24], mylab[24] , mzlab[24] ;
00251 
00252       sprintf(xlab,"%4.0f mm [%c]",-xx,(xx<0.0)?'R':'L') ;
00253       sprintf(ylab,"%4.0f mm [%c]",-yy,(yy<0.0)?'A':'P') ;
00254       sprintf(zlab,"%4.0f mm [%c]", zz,(zz<0.0)?'I':'S') ;
00255 
00256       LOAD_FVEC3(tv,xx,yy,zz);
00257       mv = THD_tta_to_mni(tv); UNLOAD_FVEC3(mv,mx,my,mz);
00258       sprintf(mxlab,"%4.0f mm [%c]",mx,(mx>=0.0)?'R':'L') ;
00259       sprintf(mylab,"%4.0f mm [%c]",my,(my>=0.0)?'A':'P') ;
00260       sprintf(mzlab,"%4.0f mm [%c]",mz,(mz< 0.0)?'I':'S') ;
00261 
00262       rbuf = AFMALL(char, 500) ;
00263       sprintf(rbuf,"%s\n"
00264                    "Focus point=%s,%s,%s {T-T Atlas}\n"
00265                    "           =%s,%s,%s {MNI Brain}\n"
00266                    "\n"
00267                    "***** Not near any region stored in database *****\n" ,
00268               WAMI_HEAD , xlab,ylab,zlab , mxlab,mylab,mzlab ) ;
00269       RETURN(rbuf) ;
00270    }
00271 
00272    /*-- bubble-sort what we found, by radius --*/
00273 
00274    if( nfind > 1 ){  /* don't have to sort only 1 result */
00275      int swap, tmp ;
00276      do{
00277         swap=0 ;
00278         for( ii=1 ; ii < nfind ; ii++ ){
00279            if( rr_find[ii-1] > rr_find[ii] ){
00280              tmp = rr_find[ii-1]; rr_find[ii-1] = rr_find[ii]; rr_find[ii] = tmp;
00281              tmp = b2_find[ii-1]; b2_find[ii-1] = b2_find[ii]; b2_find[ii] = tmp;
00282              tmp = b4_find[ii-1]; b4_find[ii-1] = b4_find[ii]; b4_find[ii] = tmp;
00283              swap++ ;
00284            }
00285         }
00286      } while(swap) ;
00287    }
00288 
00289    /*-- find anatomical label for each found marker, make result string --*/
00290 
00291    INIT_SARR(sar) ; ADDTO_SARR(sar,WAMI_HEAD) ;
00292 
00293    /* 04 Apr 2002: print coordinates (LPI) as well (the HH-PB addition) */
00294 
00295    { char lbuf[128], xlab[24], ylab[24] , zlab[24] ;
00296      sprintf(xlab,"%4.0f mm [%c]",-xx,(xx<0.0)?'R':'L') ;
00297      sprintf(ylab,"%4.0f mm [%c]",-yy,(yy<0.0)?'A':'P') ;
00298      sprintf(zlab,"%4.0f mm [%c]", zz,(zz<0.0)?'I':'S') ;
00299      sprintf(lbuf,"Focus point=%s,%s,%s {T-T Atlas}",xlab,ylab,zlab) ;
00300      ADDTO_SARR(sar,lbuf) ;
00301    }
00302 
00303    /* 29 Apr 2002: print MNI coords as well */
00304 
00305    { THD_fvec3 tv , mv ;
00306      float mx,my,mz ;
00307      char mxlab[24], mylab[24] , mzlab[24] , lbuf[128] ;
00308      LOAD_FVEC3(tv,xx,yy,zz);
00309      mv = THD_tta_to_mni(tv); UNLOAD_FVEC3(mv,mx,my,mz);
00310      sprintf(mxlab,"%4.0f mm [%c]",mx,(mx>=0.0)?'R':'L') ;
00311      sprintf(mylab,"%4.0f mm [%c]",my,(my>=0.0)?'A':'P') ;
00312      sprintf(mzlab,"%4.0f mm [%c]",mz,(mz< 0.0)?'I':'S') ;
00313      sprintf(lbuf,"Focus point=%s,%s,%s {MNI Brain}\n",mxlab,mylab,mzlab) ;
00314      ADDTO_SARR(sar,lbuf) ;
00315    }
00316 
00317    rff = -1 ;  /* rff = radius of last found label */
00318 
00319    for( ff=0 ; ff < nfind ; ff++ ){
00320       b2f = b2_find[ff] ; b4f = b4_find[ff] ; b2lab = NULL ; b4lab = NULL ;
00321 
00322       if( b2f != 0 ){                               /* find label     */
00323          for( ii=0 ; ii < TTO_COUNT ; ii++ )        /* in AFNI's list */
00324             if( b2f == TTO_list[ii].tdval ) break ;
00325          if( ii < TTO_COUNT )                       /* always true? */
00326             b2lab = TTO_list[ii].name ;
00327 
00328          if( b2lab != NULL && xx < 0 && strstr(b2lab,"Left") != NULL ) /* maybe is Right */
00329             b2lab = TTO_list[ii+1].name ;
00330       }
00331 
00332       if( b4f != 0 ){
00333          for( ii=0 ; ii < TTO_COUNT ; ii++ )
00334             if( b4f == TTO_list[ii].tdval ) break ;
00335          if( ii < TTO_COUNT )
00336             b4lab = TTO_list[ii].name ;
00337          if( b4lab != NULL && xx < 0 && strstr(b4lab,"Left") != NULL )
00338             b4lab = TTO_list[ii+1].name ;
00339       }
00340 
00341       if( b2lab == NULL && b4lab == NULL ) continue ;  /* no labels? */
00342 
00343       /* make output label into lbuf */
00344 
00345       lbuf[0] = '\0' ;
00346       if( b2lab != NULL ){
00347          if( rr_find[ff] != rff ){
00348             if( rr_find[ff] > 0 )
00349               sprintf( lbuf , "Within %d mm: %s" , rr_find[ff] , b2lab ) ;
00350             else
00351               sprintf( lbuf , "Focus point: %s" , b2lab ) ;
00352          } else {
00353             sprintf( lbuf , "             %s" , b2lab ) ;
00354          }
00355 
00356          for( kk=strlen(lbuf)-1 ; kk > 0 && lbuf[kk] == '.' ; kk-- )
00357             lbuf[kk] = '\0' ;                  /* trim trailing .'s */
00358       }
00359 
00360       if( b4lab != NULL ){
00361          kk = strlen(lbuf) ;
00362          if( kk > 0 ){
00363             sprintf( lbuf+kk , " -AND- %s" , b4lab ) ;
00364          } else if( rr_find[ff] != rff ){
00365             if( rr_find[ff] > 0 )
00366               sprintf( lbuf , "Within %d mm: %s" , rr_find[ff] , b4lab ) ;
00367             else
00368               sprintf( lbuf , "Focus point: %s" , b4lab ) ;
00369          } else {
00370             sprintf( lbuf , "             %s" , b4lab ) ;
00371          }
00372 
00373          for( kk=strlen(lbuf)-1 ; kk > 0 && lbuf[kk] == '.' ; kk-- )
00374             lbuf[kk] = '\0' ;
00375       }
00376 
00377       ADDTO_SARR(sar,lbuf) ;  /* make a list of labels */
00378 
00379       rff = rr_find[ff] ;  /* save for next time around */
00380    }
00381 
00382    /*- if didn't make any label, must produce something -*/
00383 
00384    if( sar->num == 1 ){    /* shouldn't ever happen */
00385       sprintf(lbuf,"Found %d marked but unlabeled regions???\n",nfind) ;
00386       ADDTO_SARR(sar,lbuf) ;
00387    } else {
00388       ADDTO_SARR(sar,WAMI_TAIL) ;  /* cautionary tail */
00389    }
00390 
00391    /*- convert list of labels into one big multi-line string -*/
00392 
00393    for( nfind=ii=0 ; ii < sar->num ; ii++ ) nfind += strlen(sar->ar[ii]) ;
00394    rbuf = AFMALL(char, nfind + 2*sar->num + 32 ) ; rbuf[0] = '\0' ;
00395    for( ii=0 ; ii < sar->num ; ii++ ){
00396       strcat(rbuf,sar->ar[ii]) ; strcat(rbuf,"\n") ;
00397    }
00398 
00399    DESTROY_SARR(sar) ; RETURN(rbuf) ;
00400 }