1 /*** Whereami.c modified 1/11/05 -- main function by Mike Angstadt of U Chicago ***/
2
3 #define MAIN
4 #define SUMA_noFunc
5
6 #include "mrilib.h"
7 #include <stdio.h>
8 #include <stdlib.h>
9 #include "matrix.h"
10 #include "suma_objs.h" /* 21 Apr 2020 */
11 /*------------------------------------------------------------*/
12
13 /**Original code by Mike Angstadt *******************************************
14 Main function added by Mike Angstadt on 1/12/05
15 usage: whereami x y z [output format]
16 where x,y,z are float coordinates in tlrc space
17 and output format is a 0 or 1
18 0 (default) just outputs the string as it would appear from the interactive
19 AFNI Where am I? command
20 1 outputs the string as a tab-delimited list of the form:
21 Focus point: Some area <tab> Within 6 mm: Some area <tab> etc
22 **************************************************************************
23 New version by Ziad S. Saad Feb 06
24 */
25
26 #define zischar(ch) ( ( ((ch) >= 'A' && (ch) <= 'Z' ) || ((ch) >= 'a' && (ch) <= 'z' ) ) ? 1 : 0 )
27 #define isnakedarg(s) ( ( (s)[0] == '-' && strlen(s) > 1 && zischar((s)[1]) ) ? 0 : 1 )
28
29 int compute_overlap(char *bmsk, byte *cmask, int ncmask, int dobin,
30 int N_atlas_names, char **atlas_names, ATLAS_LIST *atlas_alist);
31 static float *make_coord_list(char *coord_file, int *nxyz, int dicom);
32
33
PrettyRef(char * ref)34 char *PrettyRef(char *ref) {
35 int i=0;
36 char *pj=NULL, *pea=NULL;
37 static char strbuf[500];
38
39 pj = strstr(ref, "-> ");
40 if (!pj || pj == ref) return(ref);
41
42 pea = pj; /* now go back until no more - are found */
43 while (pea[0] == '-') --pea;
44
45 pj = pj + 3; /* start of journal reference */
46
47 /* copy name of area */
48 i = 0;
49 while (ref<=pea) {
50 strbuf[i] = *ref; ++ref;++i;
51 }
52 strbuf[i] = '\0';
53
54 /* now add the reference */
55 snprintf(strbuf+i, 490*sizeof(char), "\n -> %s", pj);
56
57 return(strbuf);
58 }
59
print_atlas_reference(char * atname)60 int print_atlas_reference(char *atname)
61 {
62 int i = 0, N_refs=0;
63 char **refs=NULL;
64 ATLAS *atlas;
65 atlas = Atlas_With_Trimming(atname, 1, NULL);
66 if (atlas && ATL_COMMENT(atlas)) {
67 print_atlas_comment(atlas);
68 }
69 else {
70 refs = atlas_reference_string_list(atname, &N_refs);
71 while (i < N_refs) {
72 printf("%s", PrettyRef(refs[i]));
73 ++i;
74 }
75 if (refs) refs = free_names_list(refs, N_refs);
76 }
77
78 return(1);
79 }
80
81
whereami_usage(ATLAS_LIST * atlas_alist,int detail)82 void whereami_usage(ATLAS_LIST *atlas_alist, int detail)
83 {
84 /* print help message in three sections */
85 printf(
86 "Usage: whereami [x y z [output_format]] [-lpi/-spm] [-atlas ATLAS] \n"
87 " ++ Reports brain areas located at x y z mm in some template space\n"
88 " ++ according to atlases present with your AFNI installation.\n"
89 " ++ Show the contents of available atlases\n"
90 " ++ Extract ROIs for certain atlas regions using symbolic notation\n"
91 " ++ Report on the overlap of ROIs with Atlas-defined regions.\n"
92 "\n%s", detail ? "":"use -h or -help for more help detail.\n");
93 if (detail) {
94 printf (
95 "Options (all options are optional):\n"
96 "-----------------------------------\n"
97 " x y z [output_format] : Specifies the x y z coordinates of the \n"
98 " location probed. Coordinate are in mm and \n"
99 " assumed to be in RAI or DICOM format, unless\n"
100 " otherwise specified (see -lpi/-spm below)\n"
101 " In the AFNI viewer, coordinate format is\n"
102 " specified above the coordinates in the top-left\n"
103 " of the AFNI controller. Right click in that spot\n"
104 " to change between RAI/DICOM and LPI/SPM.\n"
105 " NOTE I:In the output, the coordinates are reported\n"
106 " in LPI, in keeping with the convention used\n"
107 " in most publications.\n"
108 " NOTE II:To go between LPI and RAI, simply flip the \n"
109 " sign of the X and Y coordinates.\n"
110 "\n"
111 " Output_format is an optional flag where:\n"
112 " 0 is for standard AFNI 'Where am I?' format.\n"
113 " 1 is for Tab separated list, meant to be \n"
114 " friendly for use in spreadsheets. \n"
115 " The default output flag is 0. You can use\n"
116 " options -tab/-classic instead of the 0/1 flag.\n"
117 " -coord_file XYZ.1D: Input coordinates are stored in file XYZ.1D\n"
118 " Use the '[ ]' column selectors to specify the\n"
119 " X,Y, and Z columns in XYZ.1D.\n"
120 " Say you ran the following 3dclust command:\n"
121 " 3dclust -1Dformat -1clip 0.3 5 3000 func+orig'[1]' > out.1D\n"
122 " You can run whereami on each cluster's center\n"
123 " of mass with:\n"
124 " whereami -coord_file out.1D'[1,2,3]' -tab\n"
125 " NOTE: You cannot use -coord_file AND specify x,y,z on\n"
126 " command line.\n"
127 " -linkrbrain: get report from linkRbrain from list of coordinates\n"
128 " only with -coord_file and -space or -dset_space\n"
129 " -linkr_type tasks/genes: report for correlation with tasks or genes\n"
130 " Default is tasks\n"
131 " -lpi/-spm: Input coordinates' orientation is in LPI or SPM format. \n"
132 " -rai/-dicom: Input coordinates' orientation is in RAI or DICOM format.\n"
133 " NOTE: The default format for input coordinates' orientation is set by \n"
134 " AFNI_ORIENT environment variable. If it is not set, then the default \n"
135 " is RAI/DICOM\n"
136 " -space SPC: Space of input coordinates.\n"
137 " SPC can be any template space name. Without a NIML table definition,\n"
138 " the space name is limited to MNI, MNI_ANAT or TLRC (the default).\n"
139 " -classic: Classic output format (output_format = 0).\n"
140 " -tab: Tab delimited output (output_format = 1). \n"
141 " Useful for spreadsheeting.\n"
142 " -atlas ATLAS: Use atlas ATLAS for the query.\n"
143 " You can use this option repeatedly to specify\n"
144 " more than one atlas. Default is all available atlases.\n");
145 if (detail > 1) {
146 printf(" ATLAS is one of:\n");
147 print_atlas_table(atlas_alist);
148 } else {
149 printf(" Use whereami -help to see all available atlases.\n");
150 }
151
152 /* third section for usage help*/
153 printf(
154 " -dset: Determine the template space to use from this reference dataset\n"
155 " Space for human data is usually TLRC, MNI, MNI_ANAT.\n"
156 " If the space is known and a reference atlas can be found, the\n"
157 " regions will be based on the coordinates from this template space.\n"
158 " -atlas_sort: Sort results by atlas (default)\n"
159 " -zone_sort | -radius_sort: Sort by radius of search\n"
160 " -old : Run whereami in the olde (Pre Feb. 06) way.\n"
161 " -show_atlas_code: Shows integer code to area label map of the atlases\n"
162 " in use. The output is not too pretty because\n"
163 " the option is for debugging use.\n"
164 " -show_atlas_region REGION_CODE: You can now use symbolic notation to\n"
165 " select atlas regions. REGION_CODE has \n"
166 " three colon-separated elements forming it:\n"
167 " Atlas_Name:Side:Area.\n"
168 " Atlas_Name: one of the atlas names listed above.\n"
169 " If you do not have a particular atlas in your AFNI\n"
170 " installation, you'll need to download it (see below).\n"
171 " Side : Either left, right or nothing(::) for bilateral.\n"
172 " Area : A string identifying an area. The string cannot contain\n"
173 " blanks. Replace blanks by '_' for example Cerebellar Vermis\n"
174 " is Cerebellar_Vermis. You can also use the abbreviated \n"
175 " version cereb_ver and the program will try to guess at \n"
176 " what you want and offer suggestions if it can't find the\n"
177 " area or if there is ambiguity. Abbreviations are formed\n"
178 " by truncating the components (chunks) of an area's name \n"
179 " (label). For example:\n"
180 " 1- TT_Daemon::ant_cing specifies the bilateral\n"
181 " anterior cingulate in the TT_Daemon atlas.\n"
182 " 2- CA_N27_ML:left:hippo specifies the left\n"
183 " hippocampus in the CA_N27_ML atlas.\n"
184 " 3- CA_N27_MPM:right:124 specifies the right\n"
185 " ROI with integer code 124 in the CA_N27_MPM atlas\n"
186 " 4- CA_N27_ML::cereb_ver seeks the Cerebellar\n"
187 " Vermis in the CA_N27_ML atlas. However there\n"
188 " many distinct areas with this name so the program\n"
189 " will return with 'potential matches' or suggestions.\n"
190 " Use the suggestions to refine your query. For example:\n"
191 " CA_N27_ML::cereb_vermis_8\n"
192 " -mask_atlas_region REGION_CODE: Same as -show_atlas_region, plus\n"
193 " write out a mask dataset of the region.\n"
194 " -prefix PREFIX: Prefix for the output mask dataset\n"
195 " -max_areas MAX_N: Set a limit on the number of distinct areas to report.\n"
196 " This option will override the value set by the environment\n"
197 " variable AFNI_WHEREAMI_MAX_FIND, which is now set to %d\n"
198 " The variable AFNI_WHEREAMI_MAX_FIND should be set in your\n"
199 " .afnirc file.\n"
200 " -max_search_radius MAX_RAD: Set a limit on the maximum searching radius when\n"
201 " reporting results. This option will override the \n"
202 " value set by the environment variable \n"
203 " AFNI_WHEREAMI_MAX_SEARCH_RAD,\n"
204 " which is now set to %f .\n"
205 " -min_prob MIN_PROB: set minimum probability to consider in probabilistic\n"
206 " atlas output. This option will overrid the value set by the\n"
207 " environment variable AFNI_WHEREAMI_PROB_MIN (default is 1E-10)\n"
208 " NOTE: You can turn off some of the whining by setting the environment \n"
209 " variable AFNI_WHEREAMI_NO_WARN\n"
210 " -debug DEBUG: Debug flag\n"
211 " -verb VERB: Same as -debug DEBUG\n"
212 "\n"
213 "Options for determining the percent overlap of ROIs with Atlas-defined areas:\n"
214 "---------------------------------------------------------------------------\n"
215 " -bmask MASK: Report on the overlap of all non-zero voxels in MASK dataset\n"
216 " with various atlas regions. NOTE: The mask itself is not binary,\n"
217 " the masking operation results in a binary mask.\n"
218 " -omask ORDERED_MASK:Report on the overlap of each ROI formed by an integral \n"
219 " value in ORDERED_MASK. For example, if ORDERED_MASK has \n"
220 " ROIs with values 1, 2, and 3, then you'll get three \n"
221 " reports, one for each ROI value. Note that -omask and\n"
222 " -bmask are mutually exclusive.\n"
223 " -cmask MASK_COMMAND: command for masking values in BINARY_MASK, \n"
224 " or ORDERED_MASK on the fly.\n"
225 " e.g. whereami -bmask JoeROIs+tlrc \\\n"
226 " -cmask '-a JoeROIs+tlrc -expr equals(a,2)'\n"
227 " Would set to 0, all voxels in JoeROIs that are not\n"
228 " equal to 2.\n"
229 " Note that this mask should form a single sub-brick,\n"
230 " and must be at the same resolution as the bmask (binary mask) or\n"
231 " the omask (the ordered mask) datasets.\n"
232 " This option follows the style of 3dmaskdump (since the\n"
233 " code for it was, uh, borrowed from there (thanks Bob!, thanks Rick!)).\n"
234 " See '3dmaskdump -help' for more information.\n"
235 "\n"
236 "Note on the reported coordinates of the Focus Point:\n"
237 "----------------------------------------------------\n"
238 "Coordinates of the Focus Point are reported in available template spaces in\n"
239 "LPI coordinate order. The three principal spaces reported are Talairach \n"
240 " (TLRC), MNI, MNI Anatomical (MNI_ANAT).\n"
241 " The TLRC coordinates follow the convention specified by the Talairach and \n"
242 " Tournoux Atlas.\n"
243 " The MNI coordinates are derived from the TLRC ones using an approximation \n"
244 " equation.\n"
245 " The MNI Anat. coordinates are a shifted version of the MNI coordinates \n"
246 " (see Eickhoff et al. 05).\n"
247 "\n"
248 " For users who do not use the NIML table method of specifying template \n"
249 " and transformations, the MNI coordinates reported here are derived from TLRC\n"
250 " by an approximate function (the Brett transform). For transformations\n"
251 " between MNI_ANAT and TLRC coordinates, the 12 piece-wise linear transformation\n"
252 " that was used to transform the MNI_ANAT N27 brain to TLRC space is also\n"
253 " used to compute the coordinates in either direction.\n"
254 " For users who do use the NIML table method, the transformations among\n"
255 " the various Talairach, MNI and MNI_ANAT spaces may be performed a variety\n"
256 " of ways. The default method uses the Brett transform for TLRC to MNI, and\n"
257 " a simple shift for MNI to MNI_ANAT.\n"
258 "\n"
259 "How To See Atlas Data In AFNI as datasets:\n"
260 "------------------------------------------\n"
261 " If you want to view the atlases in the same session\n"
262 " that you are working with, choose one of options below.\n"
263 " For the sake of illustrations, I will assume that atlases\n"
264 " reside in directory: /user/abin/\n"
265 " 1-Load the session where atlases reside on afni's command\n"
266 " line: afni ./ /user/abin\n"
267 " 2-Set AFNI's environment variable AFNI_GLOBAL_SESSION\n"
268 " to the directory where the atlases reside.\n"
269 " You can add the following to you .afnirc file:\n"
270 " AFNI_GLOBAL_SESSION = /user/abin\n"
271 " Or, for a less permanent solution, you can set this environment\n"
272 " variable in the shell you are working in with (for csh and tcsh):\n"
273 " setenv AFNI_GLOBAL_SESSION /user/abin \n"
274 " ***********\n"
275 " BE CAREFUL: Do not use the AFNI_GLOBAL_SESSION approach\n"
276 " *********** if the data in your session is not already \n"
277 " written in +tlrc space. To be safe, you must have\n"
278 " both +tlrc.HEAD and +tlrc.BRIK for all datasets\n"
279 " in that session (directory). Otherwise, if the anat parents are\n"
280 " not properly set, you can end up applying the +tlrc transform\n"
281 " from one of the atlases instead of the proper anatomical \n"
282 " parent for that session.\n"
283 "\n"
284 " Note: You can safely ignore the:\n"
285 " ** Can't find anat parent .... \n"
286 " messages for the Atlas datasets.\n"
287 "\n", Init_Whereami_Max_Find(), Init_Whereami_Max_Rad());
288 if (detail > 1) {
289 printf(
290 "Convenient Color maps For Atlas Datasets:\n"
291 "----------------------------------------\n"
292 " Color maps (color scales) for atlas dataset should automatically be used\n"
293 " when these datasets are viewed in the overlay. To manually select a\n"
294 " a specific color scale in the AFNI GUI's overlay panel:\n"
295 " o set the color map number chooser to '**' \n"
296 " o right-click on the color map's color bar and select \n"
297 " 'Choose Colorscale'\n"
298 " o pick one of: CytoArch_ROI_256, CytoArch_ROI_256_gap, ROI_32. etc.\n"
299 " o set autorange off and set the range to the number of colors \n"
300 " in the chosen map (256, 32, etc.). \n"
301 " Color map CytoArch_ROI_256_gap was created for the proper viewing\n"
302 " of the Maximum Probability Maps of the Anatomy Toolbox.\n"
303 "\n"
304 "How To See Atlas regions overlaid in the AFNI GUI:\n"
305 "--------------------------------------------------\n"
306 " To see specific atlas regions overlaid on underlay and other overlay data,\n"
307 " 1. In Overlay control panel, check \"See TT Atlas Regions\" \n"
308 " 2. Switch view to Talairach in View Panel\n"
309 " 3. Right-click on image and select \"-Atlas colors\". In the Atlas colors\n"
310 " menu, select the colors you would like and then choose Done.\n"
311 " The images need to be redrawn to see the atlas regions, for instance,\n"
312 " by changing slices. Additional help is available in the Atlas colors\n"
313 " menu.\n"
314 " For the renderer plug-in, the underlay and overlay datasets should both\n"
315 " have Talairach view datasets actually written out to disk\n"
316 " The whereami and \"Talairach to\" functions are also available by right-\n"
317 " clicking in an image window.\n\n"
318 );
319 } else {
320 printf(" Use whereami -help to see more information on:\n"
321 " Convenient Color maps For Atlas Datasets:\n"
322 " and \n"
323 " How To See Atlas regions overlaid in the AFNI GUI:\n"
324 );
325 }
326 printf(
327 "Example 1:\n"
328 "----------\n"
329 " To find a cluster center close to the top of the brain at -12,-26, 76 (LPI),\n"
330 " whereami, assuming the coordinates are in Talairach space, would report:\n"
331 " whereami -12 -26 76 -lpi\n"
332 " ++ Input coordinates orientation set by user to LPI\n"
333 " +++++++ nearby Atlas structures +++++++\n"
334 "\n"
335 " Original input data coordinates in TLRC space\n"
336 "\n"
337 " Focus point (LPI)=\n"
338 " -12 mm [L], -26 mm [P], 76 mm [S] {TLRC}\n"
339 " -12 mm [L], -31 mm [P], 81 mm [S] {MNI}\n"
340 " -13 mm [L], -26 mm [P], 89 mm [S] {MNI_ANAT}\n"
341 "\n"
342 " Atlas CA_N27_MPM: Cytoarch. Max. Prob. Maps (N27)\n"
343 " Within 4 mm: Area 6\n"
344 " Within 7 mm: Area 4a\n"
345 "\n"
346 " Atlas CA_N27_ML: Macro Labels (N27)\n"
347 " Within 1 mm: Left Paracentral Lobule\n"
348 " Within 6 mm: Left Precentral Gyrus\n"
349 " -AND- Left Postcentral Gyrus\n"
350 "\n"
351 "Example 2:\n"
352 "----------\n"
353 " To create a mask dataset of both left and right amygdala, you can do:\n"
354 " whereami -prefix amymask -mask_atlas_region 'TT_Daemon::amygdala'\n\n"
355 "\n"
356 " Note masks based on atlas regions can be specified \"on the fly\" in \n"
357 " the same way with other afni commands as a dataset name (like 3dcalc,\n"
358 " for instance), so a mask, very often, is not needed as a separate,\n"
359 " explicit dataset on the disk.\n\n"
360 "\n"
361 "Example 3:\n"
362 "----------\n"
363 " To create a mask from a FreeSurfer 'aparc' volume parcellation:\n"
364 " (This assumes you have already run @SUMA_Make_Spec_FS, and your\n"
365 " afni distribution is recent. Otherwise update afni then run:\n"
366 " @MakeLabelTable -atlasize_labeled_dset aparc.a2009s+aseg_rank.nii\n"
367 " from the SUMA/ directory for that subject.)\n"
368 " To find the region's name, try something like:\n"
369 " whereami -atlas aparc.a2009s+aseg_rank -show_atlas_code | grep -i insula\n"
370 " Or you can try this search, assuming you screwed up the spelling:\n"
371 " whereami -atlas aparc+aseg_rank -show_atlas_code | \\\n"
372 " apsearch -word insola -stdin\n"
373 " If you really screw up the spelling try:\n"
374 " whereami -atlas aparc+aseg_rank -show_atlas_code | \\\n"
375 " sed 's/[-_]/ /g' | \\\n"
376 " apsearch -word insolent -stdin\n"
377 " Pick one area then run:\n"
378 " whereami -atlas aparc.a2009s+aseg_rank \\\n"
379 " -mask_atlas_region \\\n"
380 " aparc.a2009s+aseg_rank::ctx_rh_S_circular_insula_sup\n"
381 "\n"
382 "\n");
383
384 printf(
385 " \n---------------\n"
386 " Atlas NIML tables:\n"
387 " Atlas, templates, template spaces and transforms may all now be specified\n"
388 " in a text file that follows an XML-like format, NIML. The specifications\n"
389 " for the NIML table files will be described more fully elsewhere, but an\n"
390 " overview is presented here. By default, and soon to be included with the\n"
391 " AFNI distributions, the file AFNI_atlas_spaces.niml contains entries for\n"
392 " each of the available atlases, template spaces, templates and \n"
393 " transformations. Two other additional files may be specified and changed\n"
394 " using the environment variables, AFNI_SUPP_ATLAS and AFNI_LOCAL_ATLAS.\n"
395 " It is best to examine the provided NIML table as an example for extending\n"
396 " and modifying the various atlas definitions.\n"
397 "\n"
398 " Show atlas NIML table options:\n"
399 " -show_atlases : show all available atlases\n"
400 " -show_templates : show all available templates\n"
401 " -show_spaces : show all available template spaces\n"
402 " -show_xforms : show all available xforms\n"
403 " -show_atlas_all : show all the above\n"
404 " -show_atlas_dset : print dataset associated with each atlas\n"
405 " can be used with -atlas option above"
406 "\n"
407 " -show_available_spaces srcspace : show spaces that are available from\n"
408 " the source space\n"
409 " -show_chain srcspace destspace : show the chain of transformations\n"
410 " needed to go from one space to another\n"
411 " -calc_chain srcspace destspace : compute the chain of transformations\n"
412 " combining and inverting transformations where possible\n"
413 " examples: convert coordinates from TT_N27 to MNI or MNI anat space\n"
414 " whereami -calc_chain TT_N27 MNI -xform_xyz_quiet 10 20 30\n"
415 " whereami -calc_chain TT_N27 MNI -xform_xyz_quiet 0 0 0\n"
416 " whereami -calc_chain TT_N27 MNIA -xform_xyz_quiet 0 0 0\n"
417 " -xform_xyz : used with calc_chain, takes the x,y,z coordinates and \n"
418 " applies the combined chain of transformations to compute\n"
419 " a new x,y,z coordinate\n"
420 " -xform_xyz_quiet : Same as -xform_xyz but only ouputs the final result\n"
421 " -coord_out outfile : with -xform_xyz, -coord_file and -calc_chain, \n"
422 " specifies an output file for transformed coordinates\n"
423 " If not specified, coord_files will be transformed and printed\n"
424 " to stdout \n"
425 "Note setting the environment variable AFNI_WAMI_DEBUG will show detailed\n"
426 " progress throughout the various functions called within whereami.\n"
427 " For spaces defined using a NIML table, a Dijkstra search is used to find\n"
428 " the shortest path between spaces. Each transformation carries with it a\n"
429 " distance attribute that is used for this computation. By modifying this\n"
430 " field, the user can control which transformations are preferred.\n\n"
431 " -web_atlas_type XML/browser/struct : report results from web-based atlases\n"
432 " using XML output to screen, open a browser for output or just\n"
433 " return the name of the structure at the coordinate\n"
434 " -html : put whereami output in html format for display in a browser\n"
435
436 " \n---------------\n"
437 " More information about Atlases in AFNI can be found here:\n"
438 " https://afni.nimh.nih.gov/pub/dist/doc/htmldoc/template_atlas/framework.html\n"
439 " Class document illustrating whereami usage:\n"
440 " https://afni.nimh.nih.gov/pub/dist/edu/latest/afni11_roi/afni11_roi.pdf\n"
441 "---------------\n"
442 );
443
444 if (detail > 1) {
445 printf(
446 "Global Options (available to all AFNI/SUMA programs)\n"
447 "%s\n%s",
448 SUMA_Offset_SLines(get_help_help(),2), get_gopt_help());
449 }
450
451 printf("Thanks to Kristina Simonyan for feedback and testing.\n");
452
453 }
454 PRINT_COMPILE_DATE ;
455 return;
456 }
457
458
459 /*----------------------------------------------------------------------------*/
460
main(int argc,char ** argv)461 int main(int argc, char **argv)
462 {
463 float x, y, z, xi, yi, zi;
464 char *string, *shar = NULL;
465 int output = 0;
466 int nakedland = 0, k = 0, Show_Atlas_Code=0;
467 int iarg, dicom = 1, i, nakedarg, arglen, ixyz=0, nxyz=0;
468 AFNI_ATLAS *aa = NULL;
469 AFNI_ATLAS_REGION *aar= NULL;
470 byte OldMethod = 0;
471 char **atlas_names=NULL;
472 int N_atlas_names = 0, nbest = 0;
473 byte atlas_sort = 1, write_mask=0;
474 ATLAS_SEARCH *as=NULL;
475 char *mskpref= NULL, *bmsk = NULL;
476 byte *cmask=NULL ; int ncmask=0 ;
477 int dobin = 0, N_areas, mni;
478 char *coord_file=NULL;
479 float *coord_list = NULL, rad;
480 THD_3dim_dataset *space_dset = NULL, *atlas_dset = NULL;
481 int read_niml_atlas = 0, show_atlas = 0, show_atlas_spaces = 0;
482 int show_atlas_templates = 0, show_atlas_xforms = 0,show_atlas_dset=0;
483 int show_xform_chain = 0, calc_xform_chain=0, show_avail_space=0;
484 int show_atlas_point_lists = 0;
485 int linkrbrain_output = 0;
486 /* coordinates translated to send to linkrbrain */
487 char *linkrbrain_xml = "__temp_linkrbrain.xml";
488 /* unparsed XML response from linkrbrain */
489 char *temp_linkrbrain_results = "__temp_linkrbrain_results.xml";
490 char *linkrbrain_str = NULL;
491 char *linkrbrain_corr_test = NULL;
492 int linkr_corr_type = 0; /* linkrbrain correlation for tasks by default */
493 char *srcspace=NULL, *destspace=NULL;
494 char *coord_outfile=NULL;
495
496 ATLAS_XFORM_LIST *xfl = NULL, *cxfl = NULL;
497 float xout, yout, zout;
498 int xform_xyz = 0, xform_xyz_quiet = 0;
499 int atlas_writehard = 0, atlas_readhard = 0, alv=1, wv=1;
500 ATLAS_LIST *atlas_alist=NULL, *atlas_list=NULL, *atlas_rlist=NULL;
501 byte b1;
502 int LocalHead = wami_lh();
503 float minprob;
504
505 mainENTRY("whereami main"); machdep(); AFNI_logger("whereami",argc,argv);
506
507 b1 = 0;
508 mni = -1;
509 dobin = 0;
510 ncmask=0 ;
511 cmask=NULL ;
512 mskpref = NULL;
513 bmsk = NULL;
514 write_mask = 0;
515 dicom = -1; /* not set */
516 output = 0;
517 rad = -1.0;
518 N_areas = -1;
519 OldMethod = 0; /* Leave at 0 */
520 coord_file = NULL;
521 alv=2; wv=2;
522 xform_xyz_quiet = 0;
523 iarg = 1 ; nakedarg = 0; Show_Atlas_Code = 0; shar = NULL;
524
525 set_TT_whereami_version(alv,wv);
526 if(alv<2)
527 init_custom_atlas(); /* allow for custom atlas in old framework */
528 xi = 0.0; yi=0.0, zi=0.0;
529 set_wami_web_reqtype(WAMI_WEB_STRUCT); /* set web atlas output to simple structure */
530 set_AFNI_wami_output_mode(0); /* turn off HTML formatted output */
531
532 while( iarg < argc ){
533 arglen = strlen(argv[iarg]);
534 if(!isnakedarg(argv[iarg])) {
535 /* fprintf(stderr, "Not naked!\n"); */
536 /* Parse land */
537 nakedland = 0;
538 #ifdef USE_TRACING
539 if( strncmp(argv[iarg],"-trace",5) == 0 ){
540 DBG_trace = 1 ;
541 iarg++ ; continue ;
542 }
543 if( strncmp(argv[iarg],"-TRACE",5) == 0 ){
544 DBG_trace = 2 ;
545 iarg++ ; continue ;
546 }
547 #endif
548 for (i=1;i<arglen; ++i) argv[iarg][i] = tolower(argv[iarg][i]);
549
550 if (strcmp(argv[iarg],"-spm") == 0 || strcmp(argv[iarg],"-lpi") == 0) {
551 dicom = 0;
552 ++iarg;
553 continue;
554 }
555 if (strcmp(argv[iarg],"-ca_n27_version") == 0) {
556 fprintf(stdout,"Anatomy Toolbox Version in AFNI is:\n%s\n",
557 atlas_version_string("CA_N27_MPM"));
558 return(0);
559 }
560
561 if (strcmp(argv[iarg],"-ca_n27_reference") == 0) {
562 fprintf(stdout,"References for Anatomy Toolbox %s:\n",
563 atlas_version_string("CA_N27_MPM"));
564 print_atlas_reference("CA_N27_MPM");
565 fprintf(stdout,"\n");
566 return(0);
567 }
568
569 if ( strcmp(argv[iarg],"-rai") == 0 ||
570 strcmp(argv[iarg],"-dicom") == 0) {
571 dicom = 1;
572 ++iarg;
573 continue;
574 }
575
576 if (strcmp(argv[iarg],"-h") == 0 || strcmp(argv[iarg],"-help") == 0 ) {
577 atlas_alist = get_G_atlas_list();
578 whereami_usage(atlas_alist, strlen(argv[iarg]) > 3 ? 2:1);
579 return(0);
580 continue;
581 }
582 if (strcmp(argv[iarg],"-old") == 0 ) {
583 fprintf( stderr,
584 "** Error: This option is no longer in use\n");
585 return(1);
586 OldMethod = 1;
587 ++iarg;
588 continue;
589 }
590 if (strcmp(argv[iarg],"-space") == 0) {
591 ++iarg;
592 if (iarg >= argc) {
593 fprintf( stderr,
594 "** Error: Need parameter after -space\n"); return(1);
595 }
596 if (srcspace) {
597 fprintf( stderr,
598 "** Error: Specify space of input (or output mask) with either -space or -dset, not both\n");
599 return(1);
600 }
601
602 /* use srcspace as is on commandline */
603 srcspace = argv[iarg];
604 if ( strcmp(argv[iarg],"Paxinos_Rat_2007@Elsevier")==0 )
605 srcspace = "paxinos_rat_2007@Elsevier";
606 set_out_space(srcspace); /* make output space for mask dset */
607
608 ++iarg;
609 continue;
610 }
611
612 if (strcmp(argv[iarg],"-web_atlas_type") == 0) {
613 ++iarg;
614 if (iarg >= argc) {
615 fprintf( stderr,
616 "** Error: Need parameter after -web_atlas_type\n"); return(1);
617 }
618 if(strcmp(argv[iarg],"XML")==0) {
619 set_wami_web_reqtype(WAMI_WEB_PRINT_XML);
620 }
621 else{
622 if(strcmp(argv[iarg],"browser")==0) {
623 set_wami_web_reqtype(WAMI_WEB_BROWSER);
624 }
625 else {
626 if(strcmp(argv[iarg],"struct")==0) {
627 set_wami_web_reqtype(WAMI_WEB_STRUCT);
628 }
629 else {
630 fprintf( stderr, "** Error: option not value for web_atlas_type\n");
631 return(1);
632 }
633 }
634 }
635 ++iarg;
636 continue;
637 }
638
639 if (strcmp(argv[iarg],"-dset") == 0) {
640 ++iarg;
641 if (iarg >= argc) {
642 fprintf(stderr,"** Error: Need dset after -dset\n");
643 return(1);
644 }
645 if (srcspace) {
646 fprintf( stderr,
647 "** Error: Specify space of input (or output mask) with either -space or -dset, not both\n");
648 return(1);
649 }
650
651 if (!(space_dset = THD_open_dataset (argv[iarg]))) {
652 fprintf(stderr,"** Error: Failed to open data set %s.\n",
653 argv[iarg]);
654 return(1);
655 }
656 srcspace = THD_get_space(space_dset); /* update space if necess*/
657 set_out_space(srcspace); /* make output space for mask dset */
658 ++iarg;
659 continue;
660 }
661
662
663 if (strcmp(argv[iarg],"-zone_sort") == 0 ||
664 strcmp(argv[iarg],"-radius_sort") == 0) {
665 atlas_sort = 0;
666 ++iarg;
667 continue;
668 }
669 if (strcmp(argv[iarg],"-atlas_sort") == 0 ) {
670 atlas_sort = 1;
671 ++iarg;
672 continue;
673 }
674 if (strcmp(argv[iarg],"-alv2") == 0 ) {
675 alv = 2;
676 set_TT_whereami_version(alv,wv);
677 ++iarg;
678 continue;
679 }
680 if (strcmp(argv[iarg],"-wv2") == 0 ) {
681 wv = 2;
682 set_TT_whereami_version(alv,wv);
683 ++iarg;
684 continue;
685 }
686 if (strcmp(argv[iarg],"-classic") == 0 ) {
687 output = 0;
688 ++iarg;
689 continue;
690 }
691 if (strcmp(argv[iarg],"-tab") == 0 ) {
692 output = 1;
693 ++iarg;
694 continue;
695 }
696 if (strcmp(argv[iarg],"-atlas") == 0) {
697 ++iarg;
698 if (iarg >= argc) {
699 fprintf(stderr,
700 "** Error: Need parameter after -atlas\n"); return(1);
701 }
702 atlas_names =
703 add_to_names_list (atlas_names, &N_atlas_names, argv[iarg]);
704 ++iarg;
705 continue;
706 }
707
708 if (strcmp(argv[iarg],"-show_atlas_code") == 0) {
709 Show_Atlas_Code = 1;
710 ++iarg;
711 continue;
712 }
713
714 if (strcmp(argv[iarg],"-show_atlas_dset") == 0) {
715 show_atlas_dset = 1;
716 read_niml_atlas = 1;
717 ++iarg;
718 continue;
719 }
720
721 if (strcmp(argv[iarg],"-show_atlas_region") == 0 ||
722 strcmp(argv[iarg],"-mask_atlas_region") == 0) {
723 if (strncmp(argv[iarg],"-mask", 4) == 0) write_mask = 1;
724 ++iarg;
725 if (iarg >= argc) {
726 fprintf( stderr,
727 "** Error: Need parameter after"
728 "-show_atlas_region/-mask_atlas_region\n");
729 return(1);
730 }
731 shar = argv[iarg];
732 ++iarg;
733 continue;
734 }
735
736 if (strcmp(argv[iarg],"-coord_file") == 0) {
737 ++iarg;
738 if (iarg >= argc) {
739 fprintf(stderr,
740 "** Error: Need 1D file with coordinates after -coord_file\n");
741 return(1);
742 }
743 coord_file = argv[iarg];
744 ++iarg;
745 continue;
746 }
747
748 if (strcmp(argv[iarg],"-coord_out") == 0) {
749 ++iarg;
750 if (iarg >= argc) {
751 fprintf(stderr,
752 "** Error: Need 1D file with coordinates after -coord_out\n");
753 return(1);
754 }
755 coord_outfile = argv[iarg];
756 ++iarg;
757 continue;
758 }
759
760 if (strcmp(argv[iarg],"-linkrbrain") == 0) {
761 ++iarg;
762 linkrbrain_output = 1;
763 set_AFNI_wami_output_mode(1);
764 continue;
765 }
766
767 if (strcmp(argv[iarg],"-linkr_type") == 0) {
768 ++iarg;
769 if (iarg >= argc) {
770 fprintf(stderr,
771 "** Error: Need the word \"tasks\" or \"genes\" after linkr_type\n");
772 return(1);
773 }
774 if(strcmp(argv[iarg],"tasks")==0)
775 linkr_corr_type = 0;
776 else {
777 if(strcmp(argv[iarg],"genes")==0)
778 linkr_corr_type = 1;
779 else{
780 fprintf(stderr,
781 "** Error: Need the word \"tasks\" or \"genes\" after linkr_type\n");
782 return(1);
783 }
784 }
785 ++iarg;
786 continue;
787 }
788
789 if (strcmp(argv[iarg],"-linkrbrain_corr_test") == 0) {
790 ++iarg;
791 if (iarg >= argc) {
792 fprintf(stderr,
793 "** Error: Need filename for test XML to send"
794 " after -linkrbrain_corr_test\n");
795 return(1);
796 }
797 linkrbrain_output = 1;
798 linkrbrain_corr_test = argv[iarg];
799 ++iarg;
800 continue;
801 }
802
803
804 if (strcmp(argv[iarg],"-max_areas") == 0) {
805 ++iarg;
806 if (iarg >= argc) {
807 fprintf(stderr,"** Error: Need parameter after -max_areas\n");
808 return(1);
809 }
810 N_areas = atoi(argv[iarg]);
811 if (N_areas < 1 || N_areas > 50) {
812 fprintf(stderr,
813 "** Error: -max_areas parameter must be between 1 and 50.\n");
814 return(1);
815 }
816 Set_Whereami_Max_Find(N_areas);
817 ++iarg;
818 continue;
819 }
820
821 if (strcmp(argv[iarg],"-max_search_radius") == 0) {
822 ++iarg;
823 if (iarg >= argc) {
824 fprintf(stderr,
825 "** Error: Need parameter after -max_search_radius\n");
826 return(1);
827 }
828 rad = atof(argv[iarg]);
829 if (rad < 1.0 || rad > 9.5) {
830 fprintf(stderr,
831 "** Error: -max_search_radius parameter must "
832 "be between 1.0 and 9.5.\n");
833 return(1);
834 }
835 Set_Whereami_Max_Rad(rad);
836 ++iarg;
837 continue;
838 }
839
840 if (strcmp(argv[iarg],"-min_prob") == 0) {
841 ++iarg;
842 if (iarg >= argc) {
843 fprintf(stderr,
844 "** Error: Need parameter after -min_prob\n");
845 return(1);
846 }
847 minprob = atof(argv[iarg]);
848 if (minprob <= 0.0 || minprob > 1.0) {
849 fprintf(stderr,
850 "** Error: -min_prob parameter must be greater than 0.0 and less than or equal to 1.0\n");
851 return(1);
852 }
853 set_wami_minprob(minprob);
854 ++iarg;
855 continue;
856 }
857
858 if (strcmp(argv[iarg],"-dbg") == 0 ||
859 strcmp(argv[iarg],"-debug") == 0 ||
860 strcmp(argv[iarg],"-verb") == 0) {
861 ++iarg;
862 if (iarg >= argc) {
863 fprintf(stderr,
864 "** Error: Need parameter after -debug|-dbg|-verb\n");
865 return(1);
866 }
867 set_wami_verb(MIN_PAIR(atoi(argv[iarg]), 4));
868 ++iarg;
869 continue;
870 }
871 if (strcmp(argv[iarg],"-prefix") == 0) {
872 ++iarg;
873 if (iarg >= argc) {
874 fprintf(stderr,"** Error: Need parameter after -prefix\n");
875 return(1);
876 }
877 mskpref = argv[iarg];
878 ++iarg;
879 continue;
880 }
881
882 if (strcmp(argv[iarg],"-bmask") == 0 ||
883 strcmp(argv[iarg],"-omask") == 0 ) {
884 if (strcmp(argv[iarg],"-bmask") == 0) dobin = 1;
885 else dobin = 0;
886 ++iarg;
887 if (iarg >= argc) {
888 fprintf(stderr,"** Error: Need parameter after -bmask\n");
889 return(1);
890 }
891 if (bmsk) {
892 fprintf(stderr,
893 "** Error: -bmask and -omask are mutually exclusive.\n");
894 return(1);
895 }
896 bmsk = argv[iarg];
897
898 ++iarg;
899 continue;
900 }
901
902 if( strcmp(argv[iarg],"-cmask") == 0 ){
903 if( iarg+1 >= argc ){
904 fprintf(stderr,
905 "** Error: -cmask option requires a following argument!\n");
906 exit(1) ;
907 }
908 cmask = EDT_calcmask( argv[++iarg] , &ncmask, 0 ) ;
909 if( cmask == NULL ){
910 fprintf(stderr,"** Error: Can't compute -cmask!\n"); exit(1);
911 }
912 iarg++ ;
913 continue ;
914 }
915 if( strcmp(argv[iarg],"-readniml") == 0) {
916 ATLAS_XFORM_LIST *xfl, *cxfl;
917 ATLAS_SPACE_LIST *asl=NULL;
918 /* atlas_xform *xf; */
919
920 init_global_atlas_from_niml_files();
921 /* atlas testing */
922 xfl = get_xform_chain(asl->space+0,
923 asl->space+1);
924 print_xform_list(xfl);
925 cxfl = calc_xform_list(xfl);
926 print_xform_list(cxfl);
927 free_xform_list(xfl);
928 free_xform_list(cxfl);
929 free_global_atlas_structs();
930 exit(0);
931 }
932
933 if( strcmp(argv[iarg],"-show_available_spaces") == 0) {
934 read_niml_atlas = 1;
935 ++iarg;
936 if (iarg >= argc) {
937 fprintf(stderr,
938 "** Error: Need src space name after"
939 " -show_available_spaces\n");
940 return(1);
941 }
942 srcspace = argv[iarg];
943 show_avail_space = 1;
944 ++iarg;
945 continue;
946 }
947
948 if( strcmp(argv[iarg],"-show_chain") == 0) {
949 read_niml_atlas = 1;
950 ++iarg;
951 if (iarg >= argc) {
952 fprintf(stderr,
953 "** Error: Need src and dest spaces after -show_chain\n");
954 return(1);
955 }
956 srcspace = argv[iarg];
957 ++iarg;
958 if (iarg >= argc) {
959 fprintf(stderr,
960 "** Error: Need src and dest spaces after -show_chain\n");
961 return(1);
962 }
963 destspace = argv[iarg];
964 show_xform_chain = 1;
965 ++iarg;
966 continue;
967 }
968
969 if( strcmp(argv[iarg],"-calc_chain") == 0) {
970 read_niml_atlas = 1;
971 ++iarg;
972 if (iarg >= argc) {
973 fprintf(stderr,
974 "** Error: Need src and dest spaces after -calc_chain\n");
975 return(1);
976 }
977 srcspace = argv[iarg];
978 ++iarg;
979 if (iarg >= argc) {
980 fprintf(stderr,
981 "** Error: Need src and dest spaces after -calc_chain\n");
982 return(1);
983 }
984 destspace = argv[iarg];
985 show_xform_chain = 1;
986 calc_xform_chain = 1;
987 ++iarg;
988 continue;
989 }
990
991 if( strcmp(argv[iarg],"-xform_xyz") == 0 ||
992 strcmp(argv[iarg],"-xform_xyz_quiet") == 0){
993 if (strlen(argv[iarg]) > 12) xform_xyz_quiet = 1;
994 iarg++;
995 read_niml_atlas = 1;
996 xform_xyz = 1;
997 continue ;
998 }
999
1000 if(strcmp(argv[iarg],"-show_atlases") == 0) {
1001 iarg++;
1002 read_niml_atlas = 1;
1003 show_atlas = 1;
1004 continue ;
1005 }
1006 if(strcmp(argv[iarg],"-show_templates") == 0) {
1007 iarg++;
1008 read_niml_atlas = 1;
1009 show_atlas_templates = 1;
1010 continue ;
1011 }
1012 if(strcmp(argv[iarg],"-show_spaces") == 0) {
1013 iarg++;
1014 read_niml_atlas = 1;
1015 show_atlas_spaces = 1;
1016 continue ;
1017 }
1018 if(strcmp(argv[iarg],"-show_xforms") == 0) {
1019 iarg++;
1020 read_niml_atlas = 1;
1021 show_atlas_xforms = 1;
1022 continue ;
1023 }
1024 if(strcmp(argv[iarg],"-show_atlas_point_lists") == 0) {
1025 iarg++;
1026 read_niml_atlas = 1;
1027 show_atlas_point_lists = 1;
1028 continue ;
1029 }
1030 if(strcmp(argv[iarg],"-show_atlas_all") == 0) {
1031 iarg++;
1032 read_niml_atlas = 1;
1033 show_atlas = 1;
1034 show_atlas_spaces = 1;
1035 show_atlas_templates = 1;
1036 show_atlas_xforms = 1;
1037 show_atlas_point_lists = 1;
1038 continue ;
1039 }
1040
1041 if (strcmp(argv[iarg],"-write_hardcode_atlas") ==0){
1042 atlas_writehard = 1;
1043 ++iarg;
1044 continue;
1045 }
1046
1047 if (strcmp(argv[iarg],"-read_atlas_niml") ==0){
1048 atlas_readhard = 1;
1049 ++iarg;
1050 if (!(atlas_dset = THD_open_dataset (argv[iarg]))) {
1051 fprintf(stderr,"** Error: Failed to open data set %s.\n",
1052 argv[iarg]);
1053 return(1);
1054 }
1055 continue;
1056 }
1057 if (strcmp(argv[iarg],"-html") ==0){
1058 set_AFNI_wami_output_mode(1);
1059 iarg++;
1060 continue;
1061 }
1062
1063 { /* bad news in tennis shoes */
1064 fprintf(stderr,"** Error: bad option %s\n", argv[iarg]);
1065 suggest_best_prog_option(argv[0], argv[iarg]);
1066 return 1;
1067 }
1068
1069 } else {
1070 /* xyz format */
1071 if (nakedarg && !nakedland) {
1072 fprintf(stderr,"** Error: Keep x, y, z and output options together\n"
1073 "argument %s not appropriate here.\n", argv[iarg]);
1074 return 1;
1075 }
1076 if (nakedarg == 0) { xi = atof(argv[iarg]); nakedland = 1; }
1077 else if (nakedarg == 1) yi = atof(argv[iarg]);
1078 else if (nakedarg == 2) zi = atof(argv[iarg]);
1079 else if (nakedarg == 3) output = atoi(argv[iarg]);
1080 ++nakedarg;
1081 ++iarg;
1082 continue;
1083 }
1084
1085 fprintf(stderr,"** Error: Shouldn't be here Jim! (%s)\n", argv[iarg]);
1086 return 1;
1087 }
1088
1089 /* set verbiage */
1090 LocalHead = wami_lh();
1091
1092 /* user wants to see atlas database */
1093 if(read_niml_atlas) {
1094 init_global_atlas_from_niml_files();
1095 if(show_avail_space)
1096 report_available_spaces(srcspace);
1097 if(show_xform_chain)
1098 xfl = report_xform_chain(srcspace, destspace, !xform_xyz_quiet);
1099 if(calc_xform_chain) {
1100 cxfl = calc_xform_list(xfl);
1101 if (!xform_xyz_quiet) {
1102 print_xform_list(cxfl);/* print the xforms briefly with names only */
1103 print_all_xforms(cxfl);/* print combined list transforms with data */
1104 }
1105 }
1106 if(xform_xyz) {
1107 if(!cxfl)
1108 apply_xform_chain(xfl, xi, yi, zi, &xout, &yout, &zout);
1109 else
1110 apply_xform_chain(cxfl, xi, yi, zi, &xout, &yout, &zout);
1111 /* cxfl = calc_xform_list(xfl);*/
1112
1113 if (xform_xyz_quiet) {
1114 printf("%f %f %f\n", xout,yout,zout);
1115 } else {
1116 printf("Coords in: %f, %f, %f -> Coords out: %f, %f, %f\n",
1117 xi,yi,zi,xout,yout,zout);
1118 }
1119 }
1120
1121 if (coord_file) { /* load the XYZ coordinates from a 1D file */
1122 coord_list = make_coord_list(coord_file, &nxyz, dicom);
1123 wami_xform_coords_print(coord_list, nxyz, srcspace, destspace, coord_outfile);
1124 }
1125
1126 if(xfl)
1127 free_xform_list(xfl);
1128 if(cxfl)
1129 free_xform_list(cxfl);
1130 if(show_atlas)
1131 print_atlas_table(get_G_atlas_list());
1132 /* print_atlas_list(get_G_atlas_list());*/
1133 if(show_atlas_templates)
1134 print_template_list(get_G_templates_list());
1135 if(show_atlas_spaces)
1136 print_space_list(get_G_space_list());
1137 if(show_atlas_xforms)
1138 print_all_xforms(get_G_xform_list());
1139 if(show_atlas_point_lists)
1140 print_point_lists(get_G_atlas_list());
1141 if(show_atlas_dset){
1142 if (N_atlas_names != 0) {
1143 atlas_alist = get_G_atlas_list();
1144 /* check for missing atlases and stop in case of error */
1145 for (k=0; k < N_atlas_names; ++k) {
1146 print_atlas_dset(atlas_alist,atlas_names[k]);
1147 }
1148 }
1149 }
1150 free_global_atlas_structs();
1151 exit(0);
1152 }
1153
1154
1155 /* write out all the atlases that are hard-coded in AFNI to NIML files */
1156 if(atlas_writehard) {
1157 AFNI_atlas_list_to_niml();
1158 exit(0);
1159 }
1160
1161 /* read all the atlases that are hard-coded in AFNI to NIML files */
1162 if(atlas_readhard) {
1163 ATLAS_POINT_LIST *apl;
1164 apl = dset_niml_to_atlas_list(atlas_dset);
1165 print_atlas_point_list(apl);
1166 free_atlas_point_list(apl);
1167 exit(0);
1168 }
1169
1170 if (nakedarg >= 3 && coord_file) {
1171 /* bad combo */
1172 fprintf(stderr,"** Error: Can't specify x, y, z "
1173 "coordinates on command line AND in coord_file.\n");
1174 return(1) ;
1175 }
1176
1177
1178 if (dicom == -1) {
1179 THD_coorder cord;
1180 /* try to set based on AFNI_ORIENT */
1181 THD_coorder_fill (my_getenv("AFNI_ORIENT"), &cord);
1182 if (strcmp(cord.orcode,"RAI") == 0) {
1183 if(!AFNI_wami_output_mode())
1184 fprintf(stdout,
1185 "++ Input coordinates orientation set by default rules to %s\n",
1186 cord.orcode);
1187 }else if (strcmp(cord.orcode,"LPI") == 0) {
1188 if(!AFNI_wami_output_mode())
1189 fprintf(stdout,
1190 "++ Input coordinates orientation set by default rules to %s\n",
1191 cord.orcode);
1192 }else {
1193 fprintf(stderr,"** Error: Only RAI or LPI orientations allowed\n"
1194 "default setting returned %s\n"
1195 "You need to override AFNI_ORIENT \n"
1196 "and use either -lpi or -rai\n", cord.orcode);
1197 return 1;
1198 }
1199 } else {
1200 if(!AFNI_wami_output_mode()){
1201 if (dicom == 1)
1202 fprintf(stdout,"++ Input coordinates orientation set by user to %s\n",
1203 "RAI");
1204 else if (dicom == 0)
1205 fprintf(stdout,"++ Input coordinates orientation set by user to %s\n",
1206 "LPI");
1207 else { fprintf(stderr,"** Error: Should not happen!\n"); return(1); }
1208 }
1209 }
1210
1211 atlas_alist = get_G_atlas_list(); /* get the whole atlas list */
1212 if (N_atlas_names == 0) {
1213 /* use all atlases */
1214 for (k=0; k<atlas_alist->natlases; ++k) {
1215 atlas_names = add_to_names_list(atlas_names, &N_atlas_names,
1216 Atlas_Name(&(atlas_alist->atlas[k])));
1217 }
1218 } else {
1219 /* check for missing atlases and stop in case of error */
1220 for (k=0; k < N_atlas_names; ++k) {
1221 if (!get_Atlas_Named(atlas_names[k], atlas_alist)) {
1222 ERROR_message("Atlas %s not found in list.\n", atlas_names[k]);
1223 string = suggest_Atlas_Named(atlas_names[k], atlas_alist);
1224 if (string[0] != '\0') {
1225 fprintf(stderr," Perhaps: %s is what you want?\n", string);
1226 }
1227 exit(1);
1228 }
1229 }
1230 atlas_rlist = Atlas_Names_to_List(atlas_names, N_atlas_names);
1231 if(wami_verb()){
1232 INFO_message("reduced list of atlases");
1233 print_atlas_list(atlas_rlist);
1234 }
1235 }
1236
1237 if (!N_atlas_names) {
1238 ERROR_message("Found no atlases");
1239 exit(1);
1240 }
1241
1242 if (nakedarg < 3 && !Show_Atlas_Code && !shar && !bmsk && !coord_file) {
1243 ERROR_message("Missing useful options. See full help or simple usage below");
1244 whereami_usage(atlas_alist, 0);
1245 return 1;
1246 }
1247
1248 if (LocalHead) Set_Show_Atlas_Mode(LocalHead);
1249
1250 if (Show_Atlas_Code) {
1251 for (k=0; k < N_atlas_names; ++k) {
1252 aa = Build_Atlas(atlas_names[k], atlas_alist);
1253 Show_Atlas(aa);
1254 aa = Free_Atlas(aa);
1255 }
1256 }
1257
1258 if (shar) {
1259 Set_ROI_String_Decode_Verbosity(1); /* help the user */
1260 /* Do the match business */
1261 if (!(aar = ROI_String_Decode(shar, atlas_alist))) {
1262 ERROR_message("ROI string decoding failed.");
1263 } else {
1264 if (LocalHead) {
1265 fprintf(stderr,"User seeks the following region in atlas %s:\n",
1266 aar->atlas_name);
1267 Show_Atlas_Region(aar);
1268 }
1269 /* is this an OK atlas */
1270 if (!get_Atlas_Named(aar->atlas_name, atlas_alist)) {
1271 ERROR_message("Atlas %s not found", aar->atlas_name);
1272 exit(1);
1273 }
1274 if (aar->N_chnks < 1 && aar->id <= 0) {
1275 ERROR_message("bad or empty label");
1276 exit(1);
1277 }
1278 if (!(aa = Build_Atlas(aar->atlas_name, atlas_alist))) {
1279 ERROR_message("Failed to build atlas");
1280 exit(1);
1281 }
1282
1283 if (wami_verb() > 2) Show_Atlas(aa);
1284 as = Find_Atlas_Regions(aa,aar, NULL);
1285 /* analyze the matches, remember no left/right decisions made yet,
1286 and even if labels are present,
1287 right/left sides may not have different ids in atlas... */
1288 string = Report_Found_Regions(aa, aar, as, &nbest);
1289 if (string) {
1290 fprintf(stderr,"%s\n", string);
1291 } else {
1292 ERROR_message("NULL string returned");
1293 exit(1);
1294 }
1295 /* Now we know what matches, give me a mask */
1296 if (nbest && write_mask) {
1297 int codes[3], n_codes;
1298 THD_3dim_dataset *maskset=NULL;
1299 if (nbest > 2) {
1300 ERROR_message("Should not get this");
1301 exit(1);
1302 }
1303 n_codes = 1;
1304 codes[0] = aa->reg[as->iloc[0]]->id;
1305 if (nbest == 2) {
1306 if (aa->reg[as->iloc[0]]->id != aa->reg[as->iloc[1]]->id) {
1307 n_codes = 2;
1308 codes[1] = aa->reg[as->iloc[1]]->id;
1309 }
1310 }
1311
1312 if (!(maskset = Atlas_Region_Mask(aar, codes, n_codes,
1313 atlas_alist))) {
1314 ERROR_message("Failed to create mask");
1315 exit(1);
1316 } else {
1317 if(!equivalent_space(THD_get_space(maskset))){
1318 ERROR_message("Atlas does not match space requested.");
1319 exit(1);
1320 }
1321 tross_Make_History( "whereami" , argc, argv , maskset ) ;
1322 if (mskpref) {
1323 EDIT_dset_items( maskset,
1324 ADN_prefix , mskpref,
1325 ADN_none ) ;
1326 }
1327 if( THD_deathcon() && THD_is_file(DSET_HEADNAME(maskset)) ) {
1328 ERROR_message("Output dataset %s already exists -- "
1329 "can't overwrite", DSET_HEADNAME(maskset)) ;
1330 exit(1);
1331 }
1332
1333 if (LocalHead) {
1334 fprintf(stderr,"Writing ROI mask to %s...\n",
1335 DSET_HEADNAME(maskset));
1336 }
1337 DSET_write(maskset) ;
1338 DSET_delete(maskset); maskset = NULL;
1339 }
1340 }
1341 aar = Free_Atlas_Region(aar);
1342 if (as) as = Free_Atlas_Search(as);
1343 if (string) free(string); string = NULL;
1344 }
1345 aa = Free_Atlas(aa);
1346 }
1347
1348 /* le bmask business */
1349 if (bmsk) {
1350 compute_overlap(bmsk, cmask, ncmask, dobin, N_atlas_names,
1351 atlas_names, atlas_alist);
1352 }
1353
1354 if(cmask) free(cmask); cmask = NULL; /* Should not be needed beyond here */
1355
1356
1357 if (nakedarg < 3 && !coord_file) { /* nothing left to do */
1358 return(0);
1359 }
1360
1361 if (coord_file) { /* load the XYZ coordinates from a 1D file */
1362 coord_list = make_coord_list(coord_file, &nxyz, dicom);
1363 }
1364 else {
1365 coord_list = (float *)calloc(3, sizeof(float));
1366 if (!dicom) {
1367 coord_list[0] = -xi; coord_list[1] =-yi;
1368 }
1369 else {
1370 coord_list[0] = xi; coord_list[1] = yi;
1371 }
1372 coord_list[2] = zi;
1373 nxyz = 1;
1374 }
1375
1376
1377 if (!coord_list) {
1378 fprintf(stderr,"** Error: No coords!\n");
1379 return(1) ;
1380 }
1381
1382 #if 0
1383 /* linkRbrain output */
1384 if (linkrbrain_output){
1385 /* if just testing parsing linkrbrain output, just parse and exit*/
1386 /* for debugging only really */
1387 if(linkrbrain_corr_test){
1388 linkrbrain_XML_simple_report(linkrbrain_corr_test, 0);
1389 exit(0);
1390 }
1391
1392 if(srcspace==NULL) {
1393 ERROR_message("No source space specified. Use -space or -dset_space\n");
1394 exit(1);
1395 }
1396 if(coord_list==NULL) {
1397 ERROR_message("No coordinate list given. Use -coord_file\n");
1398 exit(1);
1399 }
1400 /* make XML input coordinate file for linkrbrain (convert to MNI if needed) */
1401 if(make_linkrbrain_xml(coord_list, nxyz, srcspace, "MNI", linkrbrain_xml,
1402 linkr_corr_type)!=0) {
1403 fprintf(stderr,"** Error: could not make XML file for linkrbrain\n");
1404 exit(1);
1405 }
1406 else {
1407 int lll = send_linkrbrain_xml(linkrbrain_xml, temp_linkrbrain_results) ;
1408 if( lll != 0 ){
1409 fprintf(stderr,"** Error: could not link to %s"
1410 " -- Check web connections\n",get_linkrbrain_site());
1411 exit(1);
1412 }
1413 else {
1414 /* parse the output here. For now just say it's there */
1415 printf(
1416 "Query of coordinates to %s correlations\n\n",get_linkrbrain_site());
1417 linkrbrain_XML_simple_report(temp_linkrbrain_results,
1418 linkr_corr_type);
1419 /* linkrbrain_str = parse_linkrbrain("linkrbrain.xml");*/
1420 }
1421 }
1422 exit(0);
1423 }
1424 #endif
1425
1426 for (ixyz = 0; ixyz < nxyz; ++ixyz) {
1427 x = coord_list[3*ixyz];
1428 y = coord_list[3*ixyz+1];
1429 z = coord_list[3*ixyz+2];
1430
1431 #if 0
1432 if (!dicom) {
1433 /* go from lpi to rai */
1434 x = -x;
1435 y = -y;
1436 }
1437 #endif
1438 if (!OldMethod) {
1439 /* the new whereami */
1440 if (atlas_sort) {
1441 if (output == 1) TT_whereami_set_outmode (TAB1_WAMI_ATLAS_SORT);
1442 else TT_whereami_set_outmode (CLASSIC_WAMI_ATLAS_SORT);
1443 } else {
1444 if (output == 1) TT_whereami_set_outmode (TAB1_WAMI_ZONE_SORT);
1445 else TT_whereami_set_outmode (CLASSIC_WAMI_ZONE_SORT);
1446 }
1447
1448 set_TT_whereami_version(alv,wv);
1449
1450 if(!atlas_rlist){
1451 atlas_list = env_atlas_list();
1452 if(!atlas_list) {
1453 atlas_list = atlas_alist;
1454 }
1455 }
1456 else {
1457 atlas_list = atlas_rlist; /* use reduced list */
1458 if (wami_verb() >= 2){
1459 INFO_message("Calling tt_whereami with this reduced"
1460 " list of atlases");
1461 print_atlas_list(atlas_rlist);
1462 }
1463 }
1464 if(space_dset) {
1465 if (LocalHead) INFO_message("Calling tt_whereami with space_dset");
1466 string = TT_whereami(x,y,z,
1467 THD_get_space(space_dset), atlas_list);
1468 } else {
1469 if (!srcspace)
1470 srcspace = TT_whereami_default_spc_name();
1471 if (LocalHead) INFO_message("Calling tt_whereami with srcspace %s",
1472 srcspace);
1473 string = TT_whereami(x,y,z, srcspace, atlas_list);
1474 }
1475 if (string) fprintf(stdout,"%s\n", string);
1476 else{
1477 if(!get_wami_web_found())
1478 fprintf(stdout,"whereami NULL string out.\n");
1479 }
1480 if (string) free(string); string = NULL;
1481 }
1482 } /* ixyz */
1483
1484 if (coord_list) free(coord_list); coord_list = NULL;
1485
1486 return 0;
1487 }
1488 /*----------------------------------------------------------------------------*/
1489 /* End whereami main */
1490 /*----------------------------------------------------------------------------*/
1491
1492 int
compute_overlap(char * bmsk,byte * cmask,int ncmask,int dobin,int N_atlas_names,char ** atlas_names,ATLAS_LIST * atlas_alist)1493 compute_overlap(char *bmsk, byte *cmask, int ncmask, int dobin,
1494 int N_atlas_names, char **atlas_names, ATLAS_LIST *atlas_alist)
1495 {
1496 byte *bmask_vol = NULL, *bba = NULL;
1497 short *ba = NULL;
1498 THD_3dim_dataset *mset=NULL, *mset_orig = NULL, *rset = NULL;
1499 ATLAS *atlas=NULL;
1500 int isb, nvox_in_mask=0, *count = NULL, dset_kind;
1501 int *ics=NULL, *unq=NULL, n_unq=0, iroi=0, nonzero, i, k;
1502 float frac=0.0, sum = 0.0, *fba=NULL;
1503 char tmps[20];
1504
1505 /* load the mask dset */
1506 if (!(mset_orig = THD_open_dataset (bmsk))) {
1507 fprintf(stderr,"** Error: Failed to open mask set %s.\n", bmsk);
1508 return(1);
1509 }
1510
1511 #if 0 /* No longer enforced here. See is_identity_xform_chain below*/
1512 /* are we in TLRC land? */
1513 if (mset_orig->view_type != VIEW_TALAIRACH_TYPE) {
1514 fprintf( stderr,
1515 "** Error: Mask set %s is not of the Talairach persuasion.\n",
1516 bmsk);
1517 return(1);
1518 }
1519 #endif
1520
1521 if (cmask) {
1522 if (ncmask != DSET_NVOX(mset_orig)) {
1523 fprintf(stderr,
1524 "** Error: Voxel number mismatch between -bmask and -cmask input.\n"
1525 "Make sure both volumes have the same number of voxels.\n");
1526
1527 return(1);
1528 }
1529 }
1530
1531 if (dobin) { /* one pass, do all */
1532 fprintf(stdout,"++ In binary mode ...\n");
1533 n_unq = 1;
1534 unq = NULL;
1535 } else {
1536 fprintf(stdout,"++ In ordered mode ...\n");
1537 /* get unique values*/
1538 unq = THD_unique_vals( mset_orig , 0, &n_unq, cmask );
1539 if (unq) {
1540 fprintf(stdout,"++ Have %d unique values of:\n", n_unq );
1541 for (k=0; k<n_unq; ++k) fprintf (stdout, " %d", unq [k]);
1542 fprintf (stdout, "\n");
1543 } else {
1544 fprintf(stdout,"++ Failed to find unique values.\n");
1545 return(1);
1546 }
1547 }
1548
1549 for (iroi=0; iroi<n_unq; ++iroi) {
1550 if (dobin) {
1551 mset = mset_orig;
1552 /* turn the mask dataset to zeros and 1s */
1553 if ((nonzero = THD_makedsetmask( mset , 0 , 1.0, 0.0 , cmask)) < 0) { /* get all non-zero values */
1554 fprintf(stderr,"** Error: No mask for you.\n");
1555 return(1);
1556 }
1557 } else {
1558 if (unq[iroi] == 0) { /* skip nonesense */
1559 fprintf(stdout,"++ Skipping unique value of 0\n");
1560 continue;
1561 } else {
1562 fprintf(stdout,
1563 "++ ========================================================================\n"
1564 ) ;
1565 fprintf(stdout,"++ Processing unique value of %d\n", unq[iroi]);
1566 }
1567 mset = EDIT_full_copy(mset_orig, "tmp_ccopy");
1568 /* turn the mask dataset to zeros and 1s */
1569 if ((nonzero =
1570 THD_makedsetmask( mset , 0 , (float)unq[iroi],
1571 (float)unq[iroi] , cmask)) < 0) {
1572 /* get all non-zero values */
1573 fprintf(stderr,"** Error: No mask for you either.\n");
1574 return(1);
1575 }
1576 }
1577 fprintf(stdout,"++ %d voxels in ROI\n", nonzero);
1578
1579 /* for each atlas */
1580 for (k=0; k < N_atlas_names; ++k) {
1581 if (!(atlas = Atlas_With_Trimming(atlas_names[k], 0, atlas_alist))) {
1582 if(wami_verb()) {
1583 fprintf(stderr,"** Warning: Atlas %s could not be loaded.\n",
1584 atlas_names[k]);
1585 }
1586 continue;
1587 }
1588
1589 if (!is_identity_xform_chain(THD_get_space(mset_orig),
1590 atlas->space)) {
1591 if (wami_verb()) {
1592 fprintf(stderr,
1593 "** Error: Not ready to deal with non-Identity transform chain.\n"
1594 "Path from input in %s to atlas %s in %s is:\n" ,
1595 THD_get_space(mset_orig),
1596 Atlas_Name(atlas), atlas->space);
1597 print_xform_chain(THD_get_space(mset_orig),
1598 atlas->space);
1599 }
1600 continue;
1601 }
1602
1603 if (is_probabilistic_atlas(atlas)) {
1604 /* not appropriate, skip*/
1605 continue;
1606 }
1607 if (atlas->adh->maxkeyval < 1) {
1608 if (wami_verb()>=2)
1609 fprintf(stderr,
1610 "** Warning: Atlas %s not suitable for this application.\n",
1611 Atlas_Name(atlas));
1612 continue;
1613 }
1614
1615 /* resample mask per atlas, use linear interpolation,
1616 cut-off at 0.5 */
1617 rset = r_new_resam_dset ( mset, ATL_DSET(atlas), 0, 0, 0, NULL,
1618 MRI_LINEAR, NULL, 1, 0);
1619 if (!rset) {
1620 fprintf(stderr,"** Error: Failed to reslice!?\n"); return(1);
1621 }
1622 /* get byte mask of regions > 0.5 */
1623 if (!(bmask_vol = THD_makemask( rset , 0 , 0.5 , 2.0 ))) {
1624 /* get all non-zero values */
1625 fprintf(stderr,"** Error: No byte for you.\n");
1626 return(1);
1627 }
1628 nvox_in_mask = 0;
1629 for (i=0; i<DSET_NVOX(ATL_DSET(atlas)); ++i) {
1630 if (bmask_vol[i]) ++nvox_in_mask;
1631 }
1632 fprintf( stdout,"++ %d voxels in atlas-resampled mask\n",
1633 nvox_in_mask);
1634 /* for each sub-brick sb */
1635 for (isb=0; isb< DSET_NVALS(ATL_DSET(atlas)); ++isb) {
1636 dset_kind = DSET_BRICK_TYPE(ATL_DSET(atlas),isb);
1637 if(dset_kind == MRI_short) {
1638 ba = DSET_BRICK_ARRAY(ATL_DSET(atlas),isb); /* short type */
1639 if (!ba) {
1640 ERROR_message("Unexpected NULL array");
1641 free(bmask_vol); bmask_vol = NULL;
1642 continue;
1643 }
1644 /* Create count array for range of integral values in atlas */
1645 count = (int *)calloc(atlas->adh->maxkeyval+1, sizeof(int));
1646 for (i=0; i<DSET_NVOX(ATL_DSET(atlas)); ++i) {
1647 if (bmask_vol[i] &&
1648 ba[i] >= atlas->adh->minkeyval) ++count[ba[i]];
1649 }
1650 }
1651 else if(dset_kind == MRI_byte) {
1652 bba = DSET_BRICK_ARRAY(ATL_DSET(atlas),isb); /* byte array */
1653 if (!bba) {
1654 ERROR_message("Unexpected NULL array");
1655 free(bmask_vol); bmask_vol = NULL;
1656 continue;
1657 }
1658 /* Create count array for range of integral values in atlas */
1659 count = (int *)calloc(atlas->adh->maxkeyval+1, sizeof(int));
1660 for (i=0; i<DSET_NVOX(ATL_DSET(atlas)); ++i) {
1661 if (bmask_vol[i] &&
1662 bba[i] >= atlas->adh->minkeyval) ++count[bba[i]];
1663 }
1664 }
1665 else if(dset_kind == MRI_float) {
1666 fba = DSET_BRICK_ARRAY(ATL_DSET(atlas),isb); /* float array */
1667 if (!fba) {
1668 ERROR_message("Unexpected NULL array");
1669 free(bmask_vol); bmask_vol = NULL;
1670 continue;
1671 }
1672 /* Create count array for range of integral values in atlas */
1673 count = (int *)calloc(atlas->adh->maxkeyval+1, sizeof(int));
1674 for (i=0; i<DSET_NVOX(ATL_DSET(atlas)); ++i) {
1675 if (bmask_vol[i] &&
1676 fba[i] >= atlas->adh->minkeyval) ++count[(int)fba[i]];
1677 }
1678 }
1679 else {
1680 ERROR_message(
1681 "Atlas %s is not of type short, byte, or float.\n",
1682 Atlas_Name(atlas));
1683 continue;
1684 }
1685
1686 /* Now form percentages */
1687 if (!unq) {
1688 fprintf(stdout,
1689 "Intersection of ROI (all non-zero values) with atlas %s (sb%d):\n",
1690 Atlas_Name(atlas), isb);
1691 } else {
1692 fprintf(stdout,
1693 "Intersection of ROI (valued %d) with atlas %s (sb%d):\n",
1694 unq[iroi],
1695 Atlas_Name(atlas), isb);
1696 }
1697
1698 /* sort the count */
1699 if (!(ics = z_idqsort (count, (atlas->adh->maxkeyval+1) ))) {
1700 fprintf(stderr,"** Error: Failed to sort!\n");
1701 return(1);
1702 }
1703
1704 sum = 0.0;
1705 for (i=0; i<=atlas->adh->maxkeyval; ++i) {
1706 if (count[i]) {
1707 if(ics[i]==0) continue; /* don't count codes of 0 */
1708 if(strcmp(
1709 STR_PRINT(Atlas_Val_Key_to_Val_Name(atlas, ics[i])),
1710 "NULL")==0)
1711 continue; /* don't count unlabeled codes */
1712 frac = (float)count[i]/(float)nvox_in_mask;
1713 sum += frac;
1714 sprintf(tmps, "%3.1f", frac*100.0);
1715 fprintf(stdout, " %-5s%% overlap with %s, code %d\n",
1716 tmps,
1717 STR_PRINT(Atlas_Val_Key_to_Val_Name(atlas, ics[i])),
1718 ics[i] );
1719 }
1720 }
1721 sprintf(tmps, "%3.1f", sum*100.0);
1722 fprintf(stdout, " -----\n"
1723 " %-5s%% of cluster accounted for.\n"
1724 "\n", tmps);
1725 /* done with count */
1726 if (count) free(count); count = NULL;
1727 if (ics) free(ics); ics = NULL;
1728 }
1729 /* done with resampled mset */
1730 DSET_delete(rset); rset = NULL;
1731 }
1732
1733 /* delete mset if not same as mset_orig */
1734 if (mset != mset_orig) DSET_delete(mset); mset = NULL;
1735 } /* iroi */
1736
1737 /* free unique values list, nothing done */
1738 if (unq) free(unq); unq = NULL;
1739
1740 /* done with mset_orig */
1741 DSET_delete(mset_orig); mset_orig = NULL;
1742
1743 return(0);
1744 }
1745
1746 /* make coordinate list from file */
make_coord_list(char * coord_file,int * nxyz,int dicom)1747 static float *make_coord_list(char *coord_file, int *nxyz, int dicom)
1748 {
1749 MRI_IMAGE * XYZ_im=NULL;
1750 float *XYZv = NULL;
1751 float *fcoord_list = NULL;
1752 int ixyz;
1753
1754
1755 XYZ_im = mri_read_1D( coord_file ) ;
1756 if( XYZ_im == NULL ){
1757 fprintf(stderr,"** Error: Can't read XYZ.1D file %s\n",coord_file);
1758 return(NULL) ;
1759 }
1760 if (XYZ_im->ny != 3) {
1761 fprintf(stderr,"** Error: Need three columns as input.\n"
1762 " Found %d columns\n", XYZ_im->ny);
1763 return(NULL) ;
1764 }
1765 XYZv = MRI_FLOAT_PTR(XYZ_im) ;
1766 fcoord_list = (float *)calloc(3*XYZ_im->nx, sizeof(float));
1767 if (!fcoord_list) {
1768 fprintf(stderr,"** Error: Failed to allocate\n");
1769 return(NULL) ;
1770 }
1771 /* copy to me own vectors */
1772 *nxyz = XYZ_im->nx;
1773 for (ixyz=0; ixyz<*nxyz; ++ixyz) {
1774 if (!dicom) {
1775 fcoord_list[3*ixyz] = -XYZv[ixyz];
1776 fcoord_list[3*ixyz+1] = -XYZv[ixyz+XYZ_im->nx];
1777 }
1778 else {
1779 fcoord_list[3*ixyz] = XYZv[ixyz];
1780 fcoord_list[3*ixyz+1] = XYZv[ixyz+XYZ_im->nx];
1781 }
1782 fcoord_list[3*ixyz+2] = XYZv[ixyz+XYZ_im->nx*2];
1783 }
1784 mri_free(XYZ_im); XYZ_im = NULL;
1785
1786 return(fcoord_list);
1787 }
1788