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