/*** Whereami.c modified 1/11/05 -- main function by Mike Angstadt of U Chicago ***/
#define MAIN
#include "mrilib.h"
#include "afni.h"
#include <stdio.h>
#include <stdlib.h>
#include "thd_ttatlas_query.h"
#include "rickr/r_new_resam_dset.h"
static int have_dseTT = -1 ;
static THD_3dim_dataset * dseTT = NULL ;
static THD_3dim_dataset * dseTT_big = NULL ; /* 01 Aug 2001 */
static int have_dseCA_EZ_MPM = -1 ;
static THD_3dim_dataset * dseCA_EZ_MPM = NULL ;
static int have_dseCA_EZ_PMaps = -1 ;
static THD_3dim_dataset * dseCA_EZ_PMaps = NULL ;
static int have_dseCA_EZ_ML = -1 ;
static THD_3dim_dataset * dseCA_EZ_ML = NULL ;
static MCW_cluster * wamiclust=NULL ;
static MCW_cluster * wamiclust_CA_EZ=NULL ;
/**Original code by Mike Angstadt *******************************************
Main function added by Mike Angstadt on 1/12/05
usage: whereami x y z [output format]
where x,y,z are float coordinates in tlrc space
and output format is a 0 or 1
0 (default) just outputs the string as it would appear from the interactive
AFNI Where am I? command
1 outputs the string as a tab-delimited list of the form:
Focus point: Some area <tab> Within 6 mm: Some area <tab> etc
**************************************************************************
New version by Ziad S. Saad Feb 06
*/
#define zischar(ch) ( ( ((ch) >= 'A' && (ch) <= 'Z' ) || ((ch) >= 'a' && (ch) <= 'z' ) ) ? 1 : 0 )
#define isnakedarg(s) ( ( (s)[0] == '-' && strlen(s) > 1 && zischar((s)[1]) ) ? 0 : 1 )
char *PrettyRef(char *ref) {
int i=0;
char *pj=NULL, *pea=NULL;
static char strbuf[500];
pj = strstr(ref, "-> ");
if (!pj || pj == ref) return(ref);
pea = pj; /* now go back until no more - are found */
while (pea[0] == '-') --pea;
pj = pj + 3; /* start of journal reference */
/* copy name of area */
i = 0;
while (ref<=pea) {
strbuf[i] = *ref; ++ref;++i;
}
strbuf[i] = '\0';
/* now add the reference */
snprintf(strbuf, 490*sizeof(char), "%s\n -> %s", strbuf, pj);
return(strbuf);
}
void whereami_usage(void)
{
int i = 0;
ENTRY("whereami_usage");
printf( "Usage: whereami [x y z [output_format]] [-lpi/-spm] [-atlas ATLAS] \n"
" ++ Reports brain areas located at x y z mm in TLRC space according \n"
" to atlases present with your AFNI installation.\n"
" ++ Show the contents of available atlases\n"
" ++ Extract ROIs for certain atlas regions using symbolic notation\n"
" ++ Report on the overlap of ROIs with Atlas-defined regions.\n"
"\n"
"Options (all options are optional):\n"
"-----------------------------------\n"
" x y z [output_format] : Specifies the x y z coordinates of the \n"
" location probed. Coordinate are in mm and \n"
" assumed to be in RAI or DICOM format, unless\n"
" otherwise specified (see -lpi/-spm below)\n"
" In the AFNI viewer, coordinate format is\n"
" specified above the coordinates in the top-left\n"
" of the AFNI controller. Right click in that spot\n"
" to change between RAI/DICOM and LPI/SPM.\n"
" NOTE I:In the output, the coordinates are reported\n"
" in LPI, in keeping with the convention used\n"
" in most publications.\n"
" NOTE II:To go between LPI and RAI, simply flip the \n"
" sign of the X and Y coordinates.\n"
"\n"
" Output_format is an optional flag where:\n"
" 0 is for standard AFNI 'Where am I?' format.\n"
" 1 is for Tab separated list, meant to be \n"
" friendly for use in spreadsheets. \n"
" The default output flag is 0. You can use\n"
" options -tab/-classic instead of the 0/1 flag.\n"
" -coord_file XYZ.1D: Input coordinates are stored in file XYZ.1D\n"
" Use the '[ ]' column selectors to specify the\n"
" X,Y, and Z columns in XYZ.1D.\n"
" Say you ran the following 3dclust command:\n"
" 3dclust -1Dformat -1clip 0.3 5 3000 func+orig'[1]' > out.1D\n"
" You can run whereami on each cluster's center\n"
" of mass with:\n"
" whereami -coord_file out.1D'[1,2,3]' -tab\n"
" NOTE: You cannot use -coord_file AND specify x,y,z on\n"
" command line.\n"
" -lpi/-spm: Input coordinates' orientation is in LPI or SPM format. \n"
" -rai/-dicom: Input coordinates' orientation is in RAI or DICOM format.\n"
" NOTE: The default format for input coordinates' orientation is set by \n"
" AFNI_ORIENT environment variable. If it is not set, then the default \n"
" is RAI/DICOM\n"
" -space SPC: Space of input coordinates.\n"
" SPC can be either MNI or TLRC which is the default.\n"
" If SPC is the MNI space, the x,y,z coordinates are transformed to\n"
" TLRC space prior to whereami query.\n"
" -classic: Classic output format (output_format = 0).\n"
" -tab: Tab delimited output (output_format = 1). \n"
" Useful for spreadsheeting.\n"
" -atlas ATLAS: Use atlas ATLAS for the query.\n"
" You can use this option repeatedly to specify\n"
" more than one atlas. Default is all available atlases.\n"
" ATLAS is one of:\n"
" %-12s: Created by tracing Talairach and Tournoux brain illustrations.\n"
" Generously contributed by Jack Lancaster and Peter Fox of RIC UTHSCSA)\n"
"\n"
" %-12s: Anatomy Toolbox's atlases, some created from cytoarchitectonic \n"
" %-12s: studies of 10 human post-mortem brains (CA_N27_MPM, CA_N27_PM). \n"
" %-12s: Generously contributed by Simon Eickhoff,\n"
" %-12s: Katrin Amunts and Karl Zilles of IME, Julich, \n"
" Germany. Please take into account the references and abide by the \n"
" warning below (provided with the Anatomy toolbox) when using these \n"
" atlases:\n",
Atlas_Code_to_Atlas_Name(AFNI_TLRC_ATLAS),
Atlas_Code_to_Atlas_Name(CA_EZ_N27_MPM_ATLAS),
Atlas_Code_to_Atlas_Name(CA_EZ_N27_ML_ATLAS),
Atlas_Code_to_Atlas_Name(CA_EZ_N27_PMAPS_ATLAS),
Atlas_Code_to_Atlas_Name(CA_EZ_N27_LR_ATLAS));
i = 0;
printf( " Anatomy Toolbox Reference and Warning:\n"
" --------------------------------------\n" );
do {
printf( " %s\n" , PrettyRef(CA_EZ_REF_STR[i]));
++i;
} while (CA_EZ_REF_STR[i][0] != '\0');
/* " [1] Auditory cortex (TE 1.0, TE 1.1, TE 1.2) : Morosan et al., Neuroimage, 2001\n"
" [2] Broca's area (BA 44, BA 45) : Amunts et al., J Comp Neurol, 1999\n"
" [3] Motor cortex (BA 4a, BA 4p, BA 6) : Geyer et al., Nature, 1996 ; S. Geyer, Springer press 2003\n"
" [4] Somatosensory cortex (BA 3a, BA 3b, BA 1 BA 2) : Geyer et al., Neuroimage, 1999 + 2000 ; Grefkes et al., Neuroimage, 2001\n"
" [5] Parietal operculum / SII (OP 1, OP 2, OP 3, OP 4) : Eickhoff et al., Cerebral Cortex, 2005a,b\n"
" [6] Amygdala (CM/LB/SF), Hippocampus (FD/CA /SUB/EC/HATA) : Amunts et al., Anat Embryol, 2005\n"
" [7] Visual cortex (BA 17, BA 18) : Amunts et al., Neuroimage, 2000\n"
" Warning:\n"
" All other areas may only be used after consultation (contact S.Eickhoff@fz-juelich.de)\n" */
printf( " \n"
" See Eickhoff et al. Neuroimage 25 (2005) for more info on:\n"
" Probability Maps (CA_N27_PM)\n"
" and Maximum Probability Maps (CA_N27_MPM)\n");
printf( " ----------------------------------------------------------\n"
"\n"
" -atlas_sort: Sort results by atlas (default)\n"
" -zone_sort | -radius_sort: Sort by radius of search\n"
" -old : Run whereami in the olde (Pre Feb. 06) way.\n"
" -show_atlas_code: Shows integer code to area label map of the atlases\n"
" in use. The output is not too pretty because\n"
" the option is for debugging use.\n"
" -show_atlas_region REGION_CODE: You can now use symbolic notation to\n"
" select atlas regions. REGION_CODE has \n"
" three colon-separated elements forming it:\n"
" Atlas_Name:Side:Area.\n"
" Atlas_Name: one of the atlas names listed above.\n"
" If you do not have a particular atlas in your AFNI\n"
" installation, you'll need to download it (see below).\n"
" Side : Either left, right or nothing(::) for bilateral.\n"
" Area : A string identifying an area. The string cannot contain\n"
" blanks. Replace blanks by '_' for example Cerebellar Vermis\n"
" is Cerebellar_Vermis. You can also use the abbreviated \n"
" version cereb_ver and the program will try to guess at \n"
" what you want and offer suggestions if it can't find the\n"
" area or if there is ambiguity. Abbreviations are formed\n"
" by truncating the components (chunks) of an area's name \n"
" (label). For example:\n"
" 1- TT_Daemon::ant_cing specifies the bilateral\n"
" anterior cingulate in the TT_Daemon atlas.\n"
" 2- CA_N27_ML:left:hippo specifies the left\n"
" hippocampus in the CA_N27_ML atlas.\n"
" 3- CA_N27_MPM:right:124 specifies the right\n"
" ROI with integer code 124 in the CA_N27_MPM atlas\n"
" 4- CA_N27_ML::cereb_ver seeks the Cerebellar\n"
" Vermis in the CA_N27_ML atlas. However there\n"
" many distinct areas with this name so the program\n"
" will return with 'potential matches' or suggestions.\n"
" Use the suggestions to refine your query. For example:\n"
" CA_N27_ML::cereb_vermis_8\n"
" -mask_atlas_region REGION_CODE: Same as -show_atlas_region, plus\n"
" write out a mask dataset of the region.\n"
" -prefix PREFIX: Prefix for the output mask dataset\n"
" -max_areas MAX_N: Set a limit on the number of distinct areas to report.\n"
" This option will override the value set by the environment\n"
" variable AFNI_WHEREAMI_MAX_FIND, which is now set to %d\n"
" The variable AFNI_WHEREAMI_MAX_FIND should be set in your\n"
" .afnirc file.\n"
" -max_search_radius MAX_RAD: Set a limit on the maximum searching radius when\n"
" reporting results. This option will override the \n"
" value set by the environment variable \n"
" AFNI_WHEREAMI_MAX_SEARCH_RAD,\n"
" which is now set to %f .\n"
" NOTE: You can turn off some of the whining by setting the environment \n"
" variable AFNI_WHEREAMI_NO_WARN\n"
" -debug DEBUG: Debug flag\n"
" -CA_N27_version: Output the version of the Anatomy Toolbox atlases and quit.\n"
" If you get warnings that AFNI's version differs from that \n"
" of the atlas' datasets then you will need to download the \n"
" latest atlas datasets from AFNI's website. You cannot use \n"
" older atlases because the atlas' integer-code to area-label\n"
" map changes from one version to the next.\n"
" To get the version of the atlas' datasets, run 3dNotes \n"
" on the atlases and look for 'Version' in one of the notes\n"
" printed out.\n"
"\n"
"Options for determining the percent overlap of ROIs with Atlas-defined areas:\n"
"---------------------------------------------------------------------------\n"
" -bmask MASK: Report on the overlap of all non-zero voxels in MASK dataset\n"
" with various atlas regions. NOTE: The mask itself is not binary,\n"
" the masking operation results in a binary mask.\n"
" -omask ORDERED_MASK:Report on the overlap of each ROI formed by an integral \n"
" value in ORDERED_MASK. For example, if ORDERED_MASK has \n"
" ROIs with values 1, 2, and 3, then you'll get three \n"
" reports, one for each ROI value. Note that -omask and\n"
" -bmask are mutually exculsive.\n"
" -cmask MASK_COMMAND: command for masking values in BINARY_MASK, \n"
" or ORDERED_MASK on the fly.\n"
" e.g. whereami -bmask JoeROIs+tlrc \\\n"
" -cmask '-a JoeROIs+tlrc -expr equals(a,2)'\n"
" Would set to 0, all voxels in JoeROIs that are not\n"
" equal to 2.\n"
" Note that this mask should form a single sub-brick,\n"
" and must be at the same resolution as BINARY_MASK or ORDERED_MASK.\n"
" This option follows the style of 3dmaskdump (since the\n"
" code for it was, uh, borrowed from there (thanks Bob!, thanks Rick!)).\n"
" See '3dmaskdump -help' for more information.\n"
"\n"
"Note on the reported coordinates of the Focus Point:\n"
"----------------------------------------------------\n"
" Coordinates of the Focus Point are reported in 3 coordinate spaces.\n"
"The 3 spaces are Talairach (TLRC), MNI, MNI Anatomical (MNI Anat.). \n"
"All three coordinates are reported in the LPI coordinate order.\n"
" The TLRC coordinates follow the convention specified by the Talairach and \n"
" Tournoux Atlas.\n"
" The MNI coordinates are derived from the TLRC ones using an approximation \n"
" equation.\n"
" The MNI Anat. coordinates are a shifted version of the MNI coordinates \n"
" (see Eickhoff et al. 05).\n"
"\n"
" However because the MNI coordinates reported here are derived from TLRC \n"
"by an approximate function it is best to derive the MNI Anat. coordinates \n"
"in a different manner. This option is possible because the MNI Anat. \n"
"coordinates are defined relative to the single-subject N27 dataset. \n"
"MNI Anat. coordinates are thus derived via the 12 piece-wise \n"
"linear transformations used to put the MNI N27 brain in TLRC space.\n"
"\n"
"Installing Atlases:\n"
"-------------------\n"
" Atlases are stored as AFNI datasets, plus perhaps an extra file or two.\n"
" These files should be placed in a location that AFNI can find. \n"
" Let us refer to this directory as ATLAS_DIR, usually it is the same as\n"
" the directory in which AFNI's binaries (such as the program afni) reside.\n"
" At a minimum, you need the TTatlas+tlrc dataset present to activate the \n"
" AFNI 'whereami' feature. To install it, if you do not have it already, \n"
" download TTatlas+tlrc* from this link: \n"
" http://afni.nimh.nih.gov/pub/dist/tgz/\n"
" and move TTatlas+tlrc* to ATLAS_DIR.\n"
" The Anatomy Toolbox atlases are in archives called CA_EZ_v*.tgz with *\n"
" indicating a particular version number. Download the archive from:\n"
" http://afni.nimh.nih.gov/pub/dist/tgz/, unpack it and move all the \n"
" files in the upacked directory into ATLAS_DIR.\n"
"\n"
"How To See Atlas Data In AFNI:\n"
"------------------------------\n"
" If you want to view the atlases in the same session\n"
" that you are working with, choose one of options below.\n"
" For the sake of illustrations, I will assume that atlases\n"
" reside in directory: /user/abin/\n"
" 1-Load the session where atlases reside on afni's command\n"
" line: afni ./ /user/abin\n"
" 2-Set AFNI's environment variable AFNI_GLOBAL_SESSION\n"
" to the directory where the atlases reside.\n"
" You can add the following to you .afnirc file:\n"
" AFNI_GLOBAL_SESSION = /user/abin\n"
" Or, for a less permanent solution, you can set this environment\n"
" variable in the shell you are working in with (for csh and tcsh):\n"
" setenv AFNI_GLOBAL_SESSION /user/abin \n"
" ***********\n"
" BE CAREFUL: Do not use the AFNI_GLOBAL_SESSION approach\n"
" *********** if the data in your session is not already \n"
" written in +tlrc space. To be safe, you must have\n"
" both +tlrc.HEAD and +tlrc.BRIK for all datasets\n"
" in that session (directory). Otherwise, if the anat parents are\n"
" not properly set, you can end up applying the +tlrc transform\n"
" from one of the atlases instead of the proper anatomical \n"
" parent for that session.\n"
"\n"
" Note: You can safely ignore the:\n"
" ** Can't find anat parent .... \n"
" messages for the Atlas datasets.\n"
"\n"
"Convenient Colormaps For Atlas Datasets:\n"
"----------------------------------------\n"
" New colormaps (colorscales) have been added\n"
" to AFNI to facilitate integral valued datasets\n"
" like ROIs and atlases. Here's what to do:\n"
" o set the color map number chooser to '**' \n"
" o right-click on the color map and select 'Choose Colorscale'\n"
" o pick one of: CytoArch_ROI_256, CytoArch_ROI_256_gap, ROI_32. etc.\n"
" o set autorange off and set the range to the number of colors \n"
" in the chosen map (256, 32, etc.). \n"
" Color maps CytoArch_ROI_256_gap was created for the proper viewing\n"
" of the Maximum Probability Maps of the Anatomy Toolbox.\n"
"\n"
"Questions Comments:\n"
"-------------------\n"
" Ziad S. Saad (ziad@nih.gov)\n"
" SSCC/NIMH/NIH/DHHS/USA\n"
"\n"
"Thanks to Kristina Simonyan for feedback and testing.\n"
"\n"
"\n",Init_Whereami_Max_Find(), Init_Whereami_Max_Rad());
EXRETURN;
}
int main(int argc, char **argv)
{
float x, y, z, xi, yi, zi, tx, ty, tz;
char *string, *fstring, atlas_name[256], *sfp=NULL, *shar = NULL;
int output = 0;
int first = 1, num = 0;
int a, nakedland = 0, k = 0, Show_Atlas_Code=0;
int iarg, dicom = 1, i, nakedarg, arglen, ixyz=0, nxyz=0;
AFNI_ATLAS *aa = NULL;
AFNI_ATLAS_REGION *aar= NULL;
int *imatch=NULL, nmatch=0;
byte isatlasused[NUMBER_OF_ATLASES];
AFNI_ATLAS_CODES ac, atlaslist[NUMBER_OF_ATLASES]={ UNKNOWN_ATLAS, UNKNOWN_ATLAS, UNKNOWN_ATLAS, UNKNOWN_ATLAS };
byte OldMethod = 0;
int N_atlaslist = 0, nbest = 0;
byte atlas_sort = 1, LocalHead = 0, write_mask=0;
ATLAS_SEARCH *as=NULL;
char *mskpref= NULL, *bmsk = NULL;
byte *cmask=NULL ; int ncmask=0 ;
int dobin = 0, N_areas, mni;
char *coord_file=NULL;
float *coord_list = NULL, rad;
THD_fvec3 tv, m;
mni = -1;
dobin = 0;
ncmask=0 ;
cmask=NULL ;
mskpref = NULL;
bmsk = NULL;
write_mask = 0;
dicom = -1; /* not set */
output = 0;
rad = -1.0;
N_areas = -1;
OldMethod = 0;
coord_file = NULL;
for (k=0; k < NUMBER_OF_ATLASES; ++k) isatlasused[k] = 0;
iarg = 1 ; nakedarg = 0; Show_Atlas_Code = 0; shar = NULL;
sprintf(atlas_name, "TT_Daemon");
xi = 0.0; yi=0.0, zi=0.0;
while( iarg < argc ){
arglen = strlen(argv[iarg]);
if(!isnakedarg(argv[iarg])) {
/* fprintf(stderr, "Not naked!\n"); */
/* Parse land */
nakedland = 0;
#ifdef USE_TRACING
if( strncmp(argv[iarg],"-trace",5) == 0 ){
DBG_trace = 1 ;
iarg++ ; continue ;
}
if( strncmp(argv[iarg],"-TRACE",5) == 0 ){
DBG_trace = 2 ;
iarg++ ; continue ;
}
#endif
for (i=1;i<arglen; ++i) argv[iarg][i] = tolower(argv[iarg][i]);
if (strcmp(argv[iarg],"-spm") == 0 || strcmp(argv[iarg],"-lpi") == 0) {
dicom = 0;
++iarg;
continue;
}
if (strcmp(argv[iarg],"-ca_n27_version") == 0) {
fprintf(stdout,"Anatomy Toolbox Version in AFNI is:\n%s\n", CA_EZ_VERSION_STR);
return(0);
}
if (strcmp(argv[iarg],"-rai") == 0 || strcmp(argv[iarg],"-dicom") == 0) {
dicom = 1;
++iarg;
continue;
}
if (strcmp(argv[iarg],"-help") == 0 ) {
whereami_usage();
return(1);
continue;
}
if (strcmp(argv[iarg],"-old") == 0 ) {
OldMethod = 1;
++iarg;
continue;
}
if (strcmp(argv[iarg],"-space") == 0) {
++iarg;
if (iarg >= argc) {
fprintf(stderr,"** Error: Need parameter after -space\n"); return(1);
}
if (strcmp(argv[iarg],"MNI") == 0 || strcmp(argv[iarg],"mni") == 0) {
mni = 1;
} else if (strcmp(argv[iarg],"TLRC") == 0 || strcmp(argv[iarg],"tlrc") == 0) {
mni = 0;
} else {
fprintf(stderr,"** Error: %s is invalid. Must use either MNI or TLRC\n", argv[iarg]);
return(1);
}
++iarg;
continue;
}
if (strcmp(argv[iarg],"-zone_sort") == 0 || strcmp(argv[iarg],"-radius_sort") == 0) {
atlas_sort = 0;
++iarg;
continue;
}
if (strcmp(argv[iarg],"-atlas_sort") == 0 ) {
atlas_sort = 1;
++iarg;
continue;
}
if (strcmp(argv[iarg],"-classic") == 0 ) {
output = 0;
++iarg;
continue;
}
if (strcmp(argv[iarg],"-tab") == 0 ) {
output = 1;
++iarg;
continue;
}
if (strcmp(argv[iarg],"-atlas") == 0) {
++iarg;
if (iarg >= argc) {
fprintf(stderr,"** Error: Need parameter after -atlas\n"); return(1);
}
if (N_atlaslist >= NUMBER_OF_ATLASES) {
fprintf(stderr,"** Error: Too many atlases specified.\n");
return(1);
}
if ( strcmp(argv[iarg],Atlas_Code_to_Atlas_Name(AFNI_TLRC_ATLAS)) == 0 ) {
if (!isatlasused[AFNI_TLRC_ATLAS]) {
isatlasused[AFNI_TLRC_ATLAS] = 1;
atlaslist[N_atlaslist] = AFNI_TLRC_ATLAS;
++N_atlaslist;
}
} else if (strcmp(argv[iarg],Atlas_Code_to_Atlas_Name(CA_EZ_N27_MPM_ATLAS) ) == 0) {
if (!isatlasused[CA_EZ_N27_MPM_ATLAS]) {
isatlasused[CA_EZ_N27_MPM_ATLAS] = 1;
atlaslist[N_atlaslist]= CA_EZ_N27_MPM_ATLAS;
++N_atlaslist;
}
} else if (strcmp(argv[iarg],Atlas_Code_to_Atlas_Name(CA_EZ_N27_ML_ATLAS)) == 0) {
if (!isatlasused[CA_EZ_N27_ML_ATLAS]) {
isatlasused[CA_EZ_N27_ML_ATLAS] = 1;
atlaslist[N_atlaslist]= CA_EZ_N27_ML_ATLAS;
++N_atlaslist;
}
} else if (strcmp(argv[iarg],Atlas_Code_to_Atlas_Name(CA_EZ_N27_PMAPS_ATLAS)) == 0) {
if (!isatlasused[CA_EZ_N27_PMAPS_ATLAS]) {
isatlasused[CA_EZ_N27_PMAPS_ATLAS] = 1;
atlaslist[N_atlaslist]= CA_EZ_N27_PMAPS_ATLAS;
++N_atlaslist;
}
} else if (strcmp(argv[iarg],Atlas_Code_to_Atlas_Name(CA_EZ_N27_LR_ATLAS)) == 0) {
if (!isatlasused[CA_EZ_N27_LR_ATLAS]) {
isatlasused[CA_EZ_N27_LR_ATLAS] = 1;
atlaslist[N_atlaslist]= CA_EZ_N27_LR_ATLAS;
++N_atlaslist;
}
} else if (strcmp(argv[iarg],"CA_EZ") == 0) {
if (!isatlasused[CA_EZ_N27_MPM_ATLAS]) {
atlaslist[N_atlaslist]= CA_EZ_N27_MPM_ATLAS;
++N_atlaslist;
}
if (!isatlasused[CA_EZ_N27_ML_ATLAS]) {
atlaslist[N_atlaslist]= CA_EZ_N27_ML_ATLAS;
++N_atlaslist;
}
if (!isatlasused[CA_EZ_N27_PMAPS_ATLAS]) {
atlaslist[N_atlaslist]= CA_EZ_N27_PMAPS_ATLAS;
++N_atlaslist;
}
isatlasused[CA_EZ_N27_MPM_ATLAS] = 1;
isatlasused[CA_EZ_N27_ML_ATLAS] = 1;
isatlasused[CA_EZ_N27_PMAPS_ATLAS] = 1;
} else {
fprintf(stderr,"** Error: Atlas name %s is not recognized\n", argv[iarg]);
return(1);
}
sprintf(atlas_name,"%s",argv[iarg]);
++iarg;
continue;
}
if (strcmp(argv[iarg],"-show_atlas_code") == 0) {
Show_Atlas_Code = 1;
++iarg;
continue;
}
if (strcmp(argv[iarg],"-show_atlas_region") == 0 || strcmp(argv[iarg],"-mask_atlas_region") == 0) {
if (strncmp(argv[iarg],"-mask", 4) == 0) write_mask = 1;
++iarg;
if (iarg >= argc) {
fprintf(stderr,"** Error: Need parameter after -show_atlas_region/-mask_atlas_region\n"); return(1);
}
shar = argv[iarg];
++iarg;
continue;
}
if (strcmp(argv[iarg],"-coord_file") == 0) {
++iarg;
if (iarg >= argc) {
fprintf(stderr,"** Error: Need parameter after -coord_file\n"); return(1);
}
coord_file = argv[iarg];
++iarg;
continue;
}
if (strcmp(argv[iarg],"-max_areas") == 0) {
++iarg;
if (iarg >= argc) {
fprintf(stderr,"** Error: Need parameter after -max_areas\n"); return(1);
}
N_areas = atoi(argv[iarg]);
if (N_areas < 1 || N_areas > 50) {
fprintf(stderr,"** Error: -max_areas parameter must be between 1 and 50.\n"); return(1);
}
Set_Whereami_Max_Find(N_areas);
++iarg;
continue;
}
if (strcmp(argv[iarg],"-max_search_radius") == 0) {
++iarg;
if (iarg >= argc) {
fprintf(stderr,"** Error: Need parameter after -max_search_radius\n"); return(1);
}
rad = atof(argv[iarg]);
if (rad < 1.0 || rad > 9.5) {
fprintf(stderr,"** Error: -max_search_radius parameter must be between 1.0 and 9.5.\n"); return(1);
}
Set_Whereami_Max_Rad(rad);
++iarg;
continue;
}
if (strcmp(argv[iarg],"-dbg") == 0 || strcmp(argv[iarg],"-debug") == 0 ) {
++iarg;
if (iarg >= argc) {
fprintf(stderr,"** Error: Need parameter after -debug\n"); return(1);
}
LocalHead = MIN_PAIR(atoi(argv[iarg]), 4);
++iarg;
continue;
}
if (strcmp(argv[iarg],"-prefix") == 0) {
++iarg;
if (iarg >= argc) {
fprintf(stderr,"** Error: Need parameter after -prefix\n"); return(1);
}
mskpref = argv[iarg];
++iarg;
continue;
}
if (strcmp(argv[iarg],"-bmask") == 0 || strcmp(argv[iarg],"-omask") == 0 ) {
if (strcmp(argv[iarg],"-bmask") == 0) dobin = 1;
else dobin = 0;
++iarg;
if (iarg >= argc) {
fprintf(stderr,"** Error: Need parameter after -bmask\n"); return(1);
}
if (bmsk) {
fprintf(stderr,"** Error: -bmask and -omask are mutually exclusive.\n"); return(1);
}
bmsk = argv[iarg];
++iarg;
continue;
}
if( strcmp(argv[iarg],"-cmask") == 0 ){
if( iarg+1 >= argc ){
fprintf(stderr,"** Error: -cmask option requires a following argument!\n");
exit(1) ;
}
cmask = EDT_calcmask( argv[++iarg] , &ncmask ) ;
if( cmask == NULL ){
fprintf(stderr,"** Error: Can't compute -cmask!\n"); exit(1);
}
iarg++ ;
continue ;
}
{ /* bad news in tennis shoes */
fprintf(stderr,"** Error: bad option %s\n", argv[iarg]);
return 1;
}
} else {
/* xyz format */
if (nakedarg && !nakedland) {
fprintf(stderr,"** Error: Keep x, y, z and output options together\n"
"argument %s not appropriate here.\n", argv[iarg]);
return 1;
}
if (nakedarg == 0) { xi = atof(argv[iarg]); nakedland = 1; }
else if (nakedarg == 1) yi = atof(argv[iarg]);
else if (nakedarg == 2) zi = atof(argv[iarg]);
else if (nakedarg == 3) output = atoi(argv[iarg]);
++nakedarg;
++iarg;
continue;
}
fprintf(stderr,"** Error: Shouldn't be here Jim! (%s)\n", argv[iarg]);
return 1;
}
if (nakedarg >= 3 && coord_file) {
/* bad combo */
fprintf(stderr,"** Error: Can't specify x, y, z coordinates on command line AND in coord_file.\n");
return(1) ;
}
if (dicom == -1) {
THD_coorder cord;
/* try to set based on AFNI_ORIENT */
THD_coorder_fill (my_getenv("AFNI_ORIENT"), &cord);
if (strcmp(cord.orcode,"RAI") == 0) {
fprintf(stdout,"++ Input coordinates orientation set by default rules to %s\n", cord.orcode);
}else if (strcmp(cord.orcode,"LPI") == 0) {
fprintf(stdout,"++ Input coordinates orientation set by default rules to %s\n", cord.orcode);
}else {
fprintf(stderr,"** Error: Only RAI or LPI orientations allowed\n"
"default setting returned %s\n"
"You need to override AFNI_ORIENT \n"
"and use either -lpi or -rai\n", cord.orcode);
return 1;
}
} else {
if (dicom == 1) fprintf(stdout,"++ Input coordinates orientation set by user to %s\n", "RAI");
else if (dicom == 0) fprintf(stdout,"++ Input coordinates orientation set by user to %s\n", "LPI");
else { fprintf(stderr,"** Error: Should not happen!\n"); return(1); }
}
if (mni == -1) {
fprintf(stdout,"++ Input coordinates space set by default rules to TLRC\n");
mni = 0;
} else if (mni == 0) {
fprintf(stdout,"++ Input coordinates space set by user to TLRC\n");
} else if (mni == 1) {
fprintf(stdout,"++ Input coordinates space set by user to MNI\n");
} else {
fprintf(stderr,"** Error: Should not happen!\n"); return(1);
}
if (N_atlaslist == 0) {
/* use all */
atlaslist[N_atlaslist] = AFNI_TLRC_ATLAS; ++N_atlaslist; isatlasused[AFNI_TLRC_ATLAS] = 1;
atlaslist[N_atlaslist] = CA_EZ_N27_MPM_ATLAS; ++N_atlaslist; isatlasused[CA_EZ_N27_MPM_ATLAS] = 1;
atlaslist[N_atlaslist] = CA_EZ_N27_ML_ATLAS; ++N_atlaslist; isatlasused[CA_EZ_N27_ML_ATLAS] = 1;
atlaslist[N_atlaslist] = CA_EZ_N27_PMAPS_ATLAS; ++N_atlaslist; isatlasused[CA_EZ_N27_PMAPS_ATLAS] = 1;
}
if (nakedarg < 3 && !Show_Atlas_Code && !shar && !bmsk && !coord_file) {
whereami_usage();
return 1;
}
if (LocalHead) Set_Show_Atlas_Mode(LocalHead);
if (Show_Atlas_Code) {
for (k=0; k < N_atlaslist; ++k) {
aa = Build_Atlas(atlaslist[k]);
Show_Atlas(aa);
aa = Free_Atlas(aa);
}
}
if (shar) {
Set_ROI_String_Decode_Verbosity(1); /* help the user */
/* Do the match business */
ac = UNKNOWN_ATLAS;
if (!(aar = ROI_String_Decode(shar, &ac))) {
ERROR_message("ROI string decoding failed.");
} else {
if (LocalHead) {
fprintf(stderr,"User seeks the following region in atlas %s:\n", Atlas_Code_to_Atlas_Name(ac));
Show_Atlas_Region(aar);
}
/* is this an OK atlas */
if (ac <= UNKNOWN_ATLAS || ac >= NUMBER_OF_ATLASES) {
ERROR_message("Atlas not found");
exit(1);
}
if (aar->N_chnks < 1 && aar->id <= 0) {
ERROR_message("bad or empty label");
exit(1);
}
if (!(aa = Build_Atlas(ac))) {
ERROR_message("Failed to build atlas");
exit(1);
}
if (LocalHead > 1) Show_Atlas(aa);
as = Find_Atlas_Regions(aa,aar, NULL);
/* analyze the matches, remember no left/right decisions made yet, and even if labels are present,
right/left sides may not have different ids in atlas... */
string = Report_Found_Regions(aa, aar, as, &nbest);
if (string) {
fprintf(stderr,"%s\n", string);
} else {
ERROR_message("NULL string returned");
exit(1);
}
/* Now we know what matches, give me a mask */
if (nbest && write_mask) {
int codes[3], n_codes;
THD_3dim_dataset *maskset=NULL;
if (nbest > 2) {
ERROR_message("Should not get this");
exit(1);
}
n_codes = 1;
codes[0] = aa->reg[as->iloc[0]]->id;
if (nbest == 2) {
if (aa->reg[as->iloc[0]]->id != aa->reg[as->iloc[1]]->id) {
n_codes = 2;
codes[1] = aa->reg[as->iloc[1]]->id;
}
}
if (!(maskset = Atlas_Region_Mask(ac, aar, codes, n_codes))) {
ERROR_message("Failed to create mask");
exit(1);
} else {
tross_Make_History( "whereami" , argc, argv , maskset ) ;
if (mskpref) {
EDIT_dset_items( maskset,
ADN_prefix , mskpref,
ADN_none ) ;
}
if( THD_deathcon() && THD_is_file(DSET_HEADNAME(maskset)) ) {
ERROR_message("Output dataset %s already exists -- can't overwrite", DSET_HEADNAME(maskset)) ;
exit(1);
}
if (LocalHead) {
fprintf(stderr,"Writing ROI mask to %s...\n", DSET_HEADNAME(maskset));
}
DSET_write(maskset) ;
DSET_delete(maskset); maskset = NULL;
}
}
aar = Free_Atlas_Region(aar);
if (as) as = Free_Atlas_Search(as);
if (string) free(string); string = NULL;
}
aa = Free_Atlas(aa);
}
/* le bmask business */
if (bmsk) {
byte *bmask_vol = NULL, *ba = NULL;
THD_3dim_dataset *mset=NULL, *mset_orig = NULL, *rset = NULL;
ATLAS_DSET_HOLDER adh;
int isb, nvox_in_mask=0, *count = NULL;
int *ics=NULL, *unq=NULL, n_unq=0, iroi=0, nonzero;
float frac=0.0, sum = 0.0;
char tmps[20];
/* load the mask dset */
if (!(mset_orig = THD_open_dataset (bmsk))) {
fprintf(stderr,"** Error: Failed to open mask set %s.\n", bmsk);
return(1);
}
/* are we in TLRC land? */
if (mset_orig->view_type != VIEW_TALAIRACH_TYPE) {
fprintf(stderr,"** Error: Mask set %s is not of the Talairach persuasion.\n", bmsk);
return(1);
}
if (cmask) {
if (ncmask != DSET_NVOX(mset_orig)) {
fprintf(stderr,"** Error: Voxel number mismatch between -bmask and -cmask input.\n"
"Make sure both volumes have the same number of voxels.\n");
return(1);
}
}
if (dobin) { /* one pass, do all */
fprintf(stdout,"++ In binary mode ...\n");
n_unq = 1;
unq = NULL;
} else {
fprintf(stdout,"++ In ordered mode ...\n");
/* get unique values*/
unq = THD_unique_vals( mset_orig , 0, &n_unq, cmask );
if (unq) {
fprintf(stdout,"++ Have %d unique values of:\n", n_unq );
for (k=0; k<n_unq; ++k) fprintf (stdout, " %d", unq [k]);
fprintf (stdout, "\n");
} else {
fprintf(stdout,"++ Failed to find unique values.\n");
return(1);
}
}
for (iroi=0; iroi<n_unq; ++iroi) {
if (dobin) {
mset = mset_orig;
/* turn the mask dataset to zeros and 1s */
if ((nonzero = THD_makedsetmask( mset , 0 , 1.0, 0.0 , cmask)) < 0) { /* get all non-zero values */
fprintf(stderr,"** Error: No mask for you.\n");
return(1);
}
} else {
if (unq[iroi] == 0) { /* skip nonesense */
fprintf(stdout,"++ Skipping unique value of 0\n");
continue;
} else {
fprintf(stdout,"++ Processing unique value of %d\n", unq[iroi]);
}
mset = EDIT_full_copy(mset_orig, "tmp_ccopy");
/* turn the mask dataset to zeros and 1s */
if ((nonzero = THD_makedsetmask( mset , 0 , (float)unq[iroi], (float)unq[iroi] , cmask)) < 0) { /* get all non-zero values */
fprintf(stderr,"** Error: No mask for you either.\n");
return(1);
}
}
fprintf(stdout,"++ %d voxels in ROI\n", nonzero);
/* for each atlas */
for (k=0; k < N_atlaslist; ++k) {
adh = Atlas_With_Trimming (atlaslist[k], 0);
if (!adh.dset) {
fprintf(stderr,"** Warning: Atlas %s could not be loaded.\n", Atlas_Code_to_Atlas_Name(atlaslist[k]));
continue;
}
if (adh.maxindexcode < 1) {
if (LocalHead) fprintf(stderr,"** Warning: Atlas %s not suitable for this application.\n", Atlas_Code_to_Atlas_Name(atlaslist[k]));
continue;
}
if (adh.maxindexcode > 255) {
fprintf(stderr,"** Warning: Max index code (%d) higher than expected.\n"
"What's cracking?.\n", adh.maxindexcode);
}
/* resample mask per atlas, use linear interpolation, cut-off at 0.5 */
rset = r_new_resam_dset ( mset, adh.dset, 0, 0, 0, NULL, MRI_LINEAR, NULL);
if (!rset) {
fprintf(stderr,"** Error: Failed to reslice!?\n"); return(1);
}
/* get byte mask of regions > 0.5 */
if (!(bmask_vol = THD_makemask( rset , 0 , 0.5 , 2.0 ))) { /* get all non-zero values */
fprintf(stderr,"** Error: No byte for you.\n");
return(1);
}
nvox_in_mask = 0;
for (i=0; i<DSET_NVOX(adh.dset); ++i) {
if (bmask_vol[i]) ++nvox_in_mask;
}
fprintf(stdout,"++ %d voxels in atlas-resampled mask\n", nvox_in_mask);
/* for each sub-brick sb */
for (isb=0; isb< DSET_NVALS(adh.dset); ++isb) {
ba = DSET_BRICK_ARRAY(adh.dset,isb);
if (!ba) { ERROR_message("Unexpected NULL array"); return(1); }
/* Create count array for range of integral values in atlas */
count = (int *)calloc(adh.maxindexcode+1, sizeof(int));
switch (adh.atcode) {
case AFNI_TLRC_ATLAS:
case CA_EZ_N27_ML_ATLAS:
case CA_EZ_N27_LR_ATLAS:
for (i=0; i<DSET_NVOX(adh.dset); ++i) {
if (bmask_vol[i] && ba[i] ) ++count[ba[i]]; /* Can't use 0 values, even if used in atlas codes */
/* such as for the AC/PC in TT_Daemon! They can't be*/
/* differentiated with this algorithm from non-brain, areas*/
}
break;
case CA_EZ_N27_MPM_ATLAS:
for (i=0; i<DSET_NVOX(adh.dset); ++i) {
if (bmask_vol[i] && ba[i] >= CA_EZ_MPM_MIN ) ++count[ba[i]];
}
break;
case CA_EZ_N27_PMAPS_ATLAS: /* not appropriate */
break;
default:
fprintf(stderr,"** Error: What is this atlas code (%d)?\n", adh.atcode);
return(1);
}
/* Now form percentages */
if (!unq) {
fprintf(stdout,"Intersection of ROI (all non-zero values) with atlas %s (sb%d):\n", Atlas_Code_to_Atlas_Name(atlaslist[k]), isb);
} else {
fprintf(stdout,"Intersection of ROI (valued %d) with atlas %s (sb%d):\n", unq[iroi], Atlas_Code_to_Atlas_Name(atlaslist[k]), isb);
}
/* sort the count */
if (!(ics = z_idqsort (count, (adh.maxindexcode+1) ))) {
fprintf(stderr,"** Error: Failed to sort!\n");
return(1);
}
sum = 0.0;
for (i=0; i<=adh.maxindexcode; ++i) {
if (count[i]) {
frac = (float)count[i]/(float)nvox_in_mask;
sum += frac;
sprintf(tmps, "%3.1f", frac*100.0);
fprintf(stdout, " %-5s%% overlap with %s, code %d\n",
tmps, STR_PRINT(Atlas_Val_to_Atlas_Name(adh, ics[i])), ics[i] );
}
}
sprintf(tmps, "%3.1f", sum*100.0);
fprintf(stdout, " -----\n"
" %-5s%% of cluster accounted for.\n"
"\n", tmps);
/* done with count */
if (count) free(count); count = NULL;
if (ics) free(ics); ics = NULL;
}
/* done with resampled mset */
DSET_delete(rset); rset = NULL;
}
/* delete mset if not same as mset_orig */
if (mset != mset_orig) DSET_delete(mset); mset = NULL;
} /* iroi */
/* free unique values list, nothing done */
if (unq) free(unq); unq = NULL;
/* done with mset_orig */
DSET_delete(mset_orig); mset_orig = NULL;
}
if(cmask) free(cmask); cmask = NULL; /* Should not be needed beyond here */
if (nakedarg < 3 && !coord_file) { /* nothing left to do */
return(0);
}
if (coord_file) { /* load the XYZ coordinates from a 1D file */
MRI_IMAGE * XYZ_im=NULL;
float *XYZv = NULL;
XYZ_im = mri_read_1D( coord_file ) ;
if( XYZ_im == NULL ){
fprintf(stderr,"** Error: Can't read XYZ.1D file %s\n",coord_file);
return(1) ;
}
if (XYZ_im->ny != 3) {
fprintf(stderr,"** Error: Need three columns as input.\n Found %d columns\n", XYZ_im->ny);
return(1) ;
}
XYZv = MRI_FLOAT_PTR(XYZ_im) ;
coord_list = (float *)calloc(3*XYZ_im->nx, sizeof(float));
if (!coord_list) {
fprintf(stderr,"** Error: Failed to allocate\n");
return(1) ;
}
/* copy to me own vectors */
nxyz = XYZ_im->nx;
for (ixyz=0; ixyz<nxyz; ++ixyz) {
coord_list[3*ixyz] = XYZv[ixyz];
coord_list[3*ixyz+1] = XYZv[ixyz+XYZ_im->nx];
coord_list[3*ixyz+2] = XYZv[ixyz+XYZ_im->nx*2];
}
mri_free(XYZ_im); XYZ_im = NULL;
} else {
coord_list = (float *)calloc(3, sizeof(float));
coord_list[0] = xi; coord_list[1] = yi; coord_list[2] = zi;
nxyz = 1;
}
if (!coord_list) {
fprintf(stderr,"** Error: No coords!\n");
return(1) ;
}
for (ixyz = 0; ixyz < nxyz; ++ixyz) {
x = coord_list[3*ixyz];
y = coord_list[3*ixyz+1];
z = coord_list[3*ixyz+2];
if (!dicom) {
/* go from lpi to rai */
x = -x;
y = -y;
}
/* coords here are now in RAI */
if (mni == 1) { /* go from mni to tlrc */
LOAD_FVEC3( tv , -x, -y, z ) ; /* next call expects input in MNI, LPI*/
m = THD_mni_to_tta( tv ); /* m units are in RAI */
if (ixyz == 0) {
fprintf(stdout,"++ Input coordinates being transformed from MNI RAI ([%.2f %.2f %.2f]) \n"
" to TLRC RAI ([%.2f %.2f %.2f]).\n",
x, y, z, m.xyz[0], m.xyz[1], m.xyz[2]);
}
x = m.xyz[0]; y = m.xyz[1]; z = m.xyz[2];
}
if (OldMethod) {
string = TT_whereami_old(x,y,z);
if (string == NULL ) { /* 30 Apr 2005 [rickr] */
fprintf(stderr,"** Error: whereami lookup failure: is TTatlas+tlrc/TTatlas.nii.gz available?\n");
fprintf(stderr," (the TTatlas+tlrc or TTatlas.nii.gz dataset must be in your PATH)\n");
return 1;
}
#if 0 /* ZSS does not work Fri Jan 20 15:54:41 EST 2006 */
if (output == 1) {
fstring = malloc(sizeof(string));
strncpy(fstring, "Focus point", 11);
num = 11;
for(a = 0; string[a] != '\0'; a++) {
/* remove header info up to Focus point:
remove newlines as you go; once you hit a *, stop */
if ((string[a] != ':') && (first == 1)) {
continue;
}
first = 0;
if ((string[a] == ' ') && (string[a-1] == ' ')) {
continue;
}
if ((string[a] == '*')) {
fstring[num] = '\0';
printf("%s\n", fstring);
break;
}
if ((string[a] != '\n')) {
if (string[a] == 'W') {
fstring[num++] = '\t';
}
fstring[num++] = string[a];
} else {
fstring[num++] = ' ';
}
}
free(fstring);
} else {
printf("%s\n", string);
}
#else
if (output == 1) {
/* ZSS: my best interpretation of the original intent of output == 1 */
fstring = strdup(string);
/* ignore everything up till Focus point */
sfp = strstr(string,"Focus point");
if (!sfp) {
fprintf(stderr,"** Error: Failed to find 'Focus point' string.\n"
"This is a beuge please inform the authors.\n");
return(1);
}
/* copy all the rest, replacing each new line followed by a non blank with a tab. */
k = 0;
while (*sfp != '\0' && sfp < string+strlen(string)) {
if (*sfp == '\n') { /* new line encountered */
/* put a tab in copy string */
fstring[k] = '\t'; ++k;
/* skip until next character */
while (!zischar(*sfp) && sfp < string+strlen(string)) ++ sfp;
}else{
/* add the character */
fstring[k] = *sfp; ++k; ++ sfp;
}
}
fstring[k] = '\0';
fprintf(stdout,"%s\n", fstring);
free(fstring); fstring = NULL;
} else {
printf("%s\n", string);
}
#endif
} /* end TT_Daemon */
if (!OldMethod) {
/* Set TT_whereami's atlas list */
for (k=0; k < N_atlaslist; ++k) {
TT_whereami_add_atlas(atlaslist[k]);
}
/* the new whereami */
if (atlas_sort) {
if (output == 1) TT_whereami_set_outmode (TAB1_WAMI_ATLAS_SORT);
else TT_whereami_set_outmode (CLASSIC_WAMI_ATLAS_SORT);
} else {
if (output == 1) TT_whereami_set_outmode (TAB1_WAMI_ZONE_SORT);
else TT_whereami_set_outmode (CLASSIC_WAMI_ZONE_SORT);
}
string = TT_whereami(x,y,z);
if (string) fprintf(stdout,"%s\n", string);
else fprintf(stdout,"whereami NULL string out.\n");
if (string) free(string); string = NULL;
}
} /* ixyz */
if (coord_list) free(coord_list); coord_list = NULL;
return 0;
}
syntax highlighted by Code2HTML, v. 0.9.1