00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068 static char g_history[] =
00069 "----------------------------------------------------------------------\n"
00070 "history:\n"
00071 "\n"
00072 "1.0 May 29, 2003 [rickr]\n"
00073 " - initial release\n"
00074 "\n"
00075 "1.1 June 11, 2003 [rickr]\n"
00076 " - small reorg of s2v_fill_mask2() (should have no effect)\n"
00077 " - improve description of -f_steps option\n"
00078 "\n"
00079 "1.2 July 21, 2003 [rickr]\n"
00080 " - make sure input points fit in output dataset\n"
00081 " - add min/max distance output, along with out-of-bounds count\n"
00082 "\n"
00083 "2.0 October 2, 2003 [rickr]\n"
00084 " - Major changes accepting surface data, surface coordinates, output data\n"
00085 " type, debug options, multiple sub-brick output, and node pair segment\n"
00086 " alterations.\n"
00087 " - Added many options: -surf_xyz_1D, -sdata_1D, -data_expr, -datum,\n"
00088 " -dnode, -dvoxel, -f_index, -f_p1_fr, -f_pn_fr,\n"
00089 " -f_p1_mm, -f_pn_mm\n"
00090 "\n"
00091 "2.1 October 20, 2003 [rickr]\n"
00092 " - call the new engine function, SUMA_LoadSpec_eng()\n"
00093 " (this will restrict the debug output from SUMA_LoadSpec())\n"
00094 "\n"
00095 "2.2 December 15, 2003 [rickr]\n"
00096 " - added program arguments '-surf_A' and '-surf_B' (-surf_A is required)\n"
00097 " - added option '-hist' (for program history)\n"
00098 " - explicit pointer init to NULL (a.o.t. memset() to 0)\n"
00099 "\n"
00100 "3.0 December 18, 2003 [rickr]\n"
00101 " - removed requirement of 2 surfaces for most functions\n"
00102 " (so now all functions work with either 1 or 2 input surfaces)\n"
00103 "\n"
00104 "3.1 January 23, 2004 [rickr]\n"
00105 " - SUMA_isINHmappable() is depricated, check with AnatCorrect field\n"
00106 " - reversed order of output for '-hist' option\n"
00107 "\n"
00108 "3.2 February 10, 2004 [rickr]\n"
00109 " - output a little more debug info for !AnatCorrect case\n"
00110 "\n"
00111 "3.3 March 26, 2004 [ziad]\n"
00112 " - DsetList added to SUMA_LoadSpec_eng() call\n"
00113 "\n"
00114 "3.4 June 21, 2004 [rickr]\n"
00115 " - Fixed -surf_xyz_1D option (broken in v3.0 on nsurf test).\n"
00116 "\n"
00117 "3.5 July 22, 2004 [rickr]\n"
00118 " - Fixed bug with test for valid sdata_1D file.\n"
00119 "\n"
00120 "3.6 July 28, 2004 [rickr]\n"
00121 " - Fixed bug: old change caused the default f_steps to revert to 1,\n"
00122 " now it is set back to 2 (thanks, Kuba).\n"
00123 "\n"
00124 "3.6a March 22, 2005 [rickr]\n"
00125 " - removed tabs\n"
00126 "----------------------------------------------------------------------\n";
00127
00128 #define VERSION "version 3.6a (March 22, 2005)"
00129
00130
00131
00132
00133
00134
00135
00136
00137 #include "mrilib.h"
00138 #include "parser.h"
00139 #include "SUMA_suma.h"
00140 #include "SUMA_3dSurf2Vol.h"
00141
00142
00143 SUMA_SurfaceViewer * SUMAg_SVv = NULL;
00144 int SUMAg_N_SVv = 0;
00145 SUMA_DO * SUMAg_DOv = NULL;
00146 int SUMAg_N_DOv = 0;
00147 SUMA_CommonFields * SUMAg_CF = NULL;
00148
00149
00150 char * gs2v_map_names[] = { "none", "mask", "mask2", "ave", "count",
00151 "min", "max", "max_abs" };
00152
00153
00154 extern void machdep( void );
00155
00156 #define MAIN
00157
00158
00159
00160 int main( int argc , char * argv[] )
00161 {
00162 SUMA_SurfSpecFile spec;
00163 node_list_t node_list;
00164 param_t params;
00165 s2v_opts_t sopt;
00166 opts_t opts;
00167 int ret_val;
00168
00169 mainENTRY("3dSurf2Vol main");
00170 machdep();
00171 AFNI_logger("3dSurf2Vol",argc,argv);
00172
00173
00174 if ( ( ret_val = init_options(&opts, argc, argv) ) != 0 )
00175 return ret_val;
00176
00177 if ( ( ret_val = validate_options(&opts, ¶ms) ) != 0 )
00178 return ret_val;
00179
00180 if ( (ret_val = set_smap_opts( &opts, ¶ms, &sopt )) != 0 )
00181 return ret_val;
00182
00183
00184 ret_val = read_surf_files(&opts, ¶ms, &spec, &sopt, &node_list);
00185
00186 if ( ret_val == 0 )
00187 ret_val = write_output( &sopt, &opts, ¶ms, &node_list, argc, argv );
00188
00189
00190 final_clean_up( &node_list );
00191
00192 return ret_val;
00193 }
00194
00195
00196
00197
00198
00199
00200 int write_output ( s2v_opts_t * sopt, opts_t * opts, param_t * p,
00201 node_list_t * N, int argc, char * argv[] )
00202 {
00203 ENTRY("write_output");
00204
00205 if ( sopt == NULL || opts == NULL || p == NULL || N == NULL )
00206 {
00207 fprintf( stderr, "** s2v_wo - bad params (%p,%p,%p,%p)\n",
00208 sopt, opts, p, N );
00209 RETURN(-1);
00210 }
00211
00212 p->oset = s2v_nodes2volume( N, p, sopt );
00213
00214 if ( p->oset == NULL )
00215 RETURN(-1);
00216
00217 EDIT_dset_items( p->oset, ADN_prefix, opts->oset_file, ADN_none );
00218
00219 if ( THD_is_file(DSET_HEADNAME(p->oset)) )
00220 {
00221 fprintf( stderr, "** cannot overwrite existing dataset '%s'\n",
00222 DSET_HEADNAME(p->oset) );
00223 DSET_delete( p->oset );
00224 RETURN(-1);
00225 }
00226
00227 tross_Copy_History( p->gpar, p->oset );
00228 tross_Make_History( PROG_NAME, argc, argv, p->oset );
00229
00230 if ( DSET_write( p->oset ) != True )
00231 {
00232 fprintf( stderr, "** failed to write dataset '%s', exiting...\n",
00233 opts->oset_file );
00234 RETURN(-1);
00235 }
00236
00237 RETURN(0);
00238 }
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254 THD_3dim_dataset * s2v_nodes2volume( node_list_t * N, param_t * p,
00255 s2v_opts_t * sopt )
00256 {
00257 THD_3dim_dataset * dout;
00258 THD_fvec3 * pary;
00259 double * ddata;
00260 void * vdata = NULL;
00261 int * cdata;
00262 float fac;
00263 int nvox, dsize, valid;
00264 int sub;
00265
00266 ENTRY("s2v_nodes2volume");
00267
00268 if ( N == NULL || ! ISVALID_DSET(p->gpar) )
00269 {
00270 fprintf(stderr,"** s2v_nodes2volume: bad params (%p,%p)\n", N, p->gpar);
00271 RETURN(NULL);
00272 }
00273
00274 dout = EDIT_empty_copy( p->gpar );
00275 if ( ! ISVALID_3DIM_DATASET( dout ) )
00276 {
00277 fprintf( stderr, "** failed EDIT_empty_copy()\n" );
00278 RETURN(NULL);
00279 }
00280
00281
00282 EDIT_dset_items( dout, ADN_nvals, p->nsubs,
00283 ADN_type, HEAD_FUNC_TYPE,
00284 ADN_func_type, FUNC_BUCK_TYPE,
00285 ADN_datum_all, sopt->datum,
00286 ADN_ntt, 0,
00287 ADN_none );
00288
00289 if ( sopt->debug > 1 )
00290 fprintf( stderr, "++ creating dataset '%s' of type '%s', nvals = %d\n",
00291 DSET_HEADNAME(dout), MRI_TYPE_name[sopt->datum], p->nsubs );
00292
00293 nvox = DSET_NVOX(p->gpar);
00294
00295
00296 if ( (ddata = (double *)malloc(nvox * sizeof(double))) == NULL )
00297 {
00298 fprintf( stderr, "** n2v: failed to allocate %d doubles\n", nvox );
00299 DSET_delete( dout );
00300 RETURN(NULL);
00301 }
00302
00303
00304 if ( (cdata = (int *)malloc(nvox * sizeof(int))) == NULL )
00305 {
00306 fprintf( stderr, "** n2v: failed to allocate %d ints\n", nvox );
00307 DSET_delete( dout ); free( ddata );
00308 RETURN(NULL);
00309 }
00310
00311
00312 if ( (pary = (THD_fvec3 *)malloc(sopt->f_steps*sizeof(THD_fvec3))) == NULL )
00313 {
00314 fprintf(stderr,"** n2v: cannot allocate %d THD_fvec3s\n",sopt->f_steps);
00315 DSET_delete( dout ); free( ddata ); free( cdata );
00316 RETURN(NULL);
00317 }
00318
00319 dsize = mri_datum_size(sopt->datum);
00320
00321 for ( sub = 0; sub < p->nsubs; sub++ )
00322 {
00323 valid = 0;
00324
00325
00326 if ( (vdata = malloc( nvox * dsize )) == NULL )
00327 {
00328 fprintf( stderr, "** failed to allocate %d bytes for vdata\n",
00329 nvox * dsize );
00330 DSET_delete( dout ); free( ddata ); free( cdata );
00331 RETURN(NULL);
00332 }
00333
00334
00335 memset(cdata, 0, nvox*sizeof(int));
00336 memset(ddata, 0, nvox*sizeof(double));
00337
00338
00339 if ( set_node_list_data( N, p, sopt, sub ) != 0 )
00340 {
00341 DSET_delete(dout);
00342 free(ddata);
00343 free(cdata);
00344 free(vdata);
00345
00346 RETURN(NULL);
00347 }
00348
00349 if ( compute_results( p, N, sopt, ddata, cdata, pary ) == 0 )
00350 valid = 1;
00351
00352 if ( ! valid )
00353 {
00354 free( ddata );
00355 free( vdata );
00356 DSET_delete( dout );
00357 RETURN(NULL);
00358 }
00359
00360
00361 fac = 0.0;
00362
00363 if ( sopt->debug > 1 )
00364 fprintf(stderr,"++ sub-brick %d, integral_doubles = %d\n",
00365 sub, integral_doubles( ddata, nvox ));
00366
00367
00368 if ( MRI_IS_INT_TYPE( sopt->datum ) && !sopt->noscale )
00369 {
00370 float amax = MCW_vol_amax( nvox, 1, 1, MRI_double, ddata );
00371
00372 if ( amax > MRI_TYPE_maxval[sopt->datum] )
00373 fac = MRI_TYPE_maxval[sopt->datum]/amax;
00374 else if ( integral_doubles( ddata, nvox ) )
00375 fac = 0.0;
00376 else if ( amax != 0.0 )
00377 fac = MRI_TYPE_maxval[sopt->datum]/amax;
00378
00379 if ( sopt->debug > 1 )
00380 fprintf(stderr,"++ fac = %f, amax = %f \n", fac, amax);
00381 }
00382
00383 EDIT_coerce_scale_type(nvox, fac, MRI_double,ddata, sopt->datum,vdata);
00384
00385
00386 EDIT_substitute_brick( dout, sub, sopt->datum, vdata );
00387 DSET_BRICK_FACTOR( dout, sub ) = (fac != 0.0) ? 1.0/fac : 0.0;
00388 }
00389
00390 if (sopt->debug > 0)
00391 fprintf(stderr,"++ %d sub-brick(s) computed\n", p->nsubs);
00392
00393 free(ddata);
00394 free(cdata);
00395 free(pary);
00396
00397 RETURN(dout);
00398 }
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408 int compute_results( param_t * p, node_list_t * N, s2v_opts_t * sopt,
00409 double * ddata, int * cdata, THD_fvec3 * pary )
00410 {
00411 THD_fvec3 p1, pn;
00412 float dist, min_dist, max_dist;
00413 int nindex, node;
00414 int oobc;
00415
00416 ENTRY("compute_results");
00417
00418 if ( !p || !N || !sopt || !ddata || !cdata )
00419 {
00420 fprintf(stderr,"** cr: bad params (%p,%p,%p,%p,%p)\n",
00421 p, N, sopt, ddata, cdata);
00422 RETURN(-1);
00423 }
00424
00425 min_dist = 9999.9;
00426 max_dist = -1.0;
00427 oobc = 0;
00428
00429 for ( nindex = 0; nindex < N->ilen; nindex++ )
00430 {
00431 node = N->ilist[nindex];
00432 if ( node < 0 || node >= N->nnodes )
00433 {
00434 fprintf(stderr,"** node %d (%d) out-of-range [0,%d]\n",
00435 nindex, node, N->nnodes);
00436 RETURN(-1);
00437 }
00438
00439
00440 p1 = N->nodes[node];
00441
00442 if ( N->depth > 1 )
00443 pn = N->nodes[node+N->nnodes];
00444 else
00445 pn = p1;
00446
00447 if ( ! sopt->sxyz_ori_gpar )
00448 {
00449 p1 = THD_dicomm_to_3dmm(p->gpar, p1);
00450 pn = THD_dicomm_to_3dmm(p->gpar, pn);
00451 }
00452
00453 if ( sopt->debug > 0 && sopt->dnode == node )
00454 {
00455 fprintf(stderr,"-- debug node: %d\n", node );
00456 fprintf(stderr,"-- orig endpts: (%f, %f, %f)\n"
00457 " (%f, %f, %f)\n",
00458 p1.xyz[0], p1.xyz[1], p1.xyz[2],
00459 pn.xyz[0], pn.xyz[1], pn.xyz[2]);
00460 }
00461
00462 adjust_endpts( sopt, &p1, &pn );
00463
00464
00465 if ( f3mm_out_of_bounds( &p1, &p->f3mm_min, &p->f3mm_max ) &&
00466 f3mm_out_of_bounds( &pn, &p->f3mm_min, &p->f3mm_max ) )
00467 {
00468 oobc++;
00469 continue;
00470 }
00471
00472 dist = dist_f3mm( &p1, &pn );
00473 if ( dist < min_dist ) min_dist = dist;
00474 if ( dist > max_dist ) max_dist = dist;
00475
00476 make_point_list( pary, &p1, &pn, sopt );
00477
00478
00479 if ( insert_list(N, p, sopt, pary, nindex, ddata, cdata) )
00480 RETURN(-1);
00481 }
00482
00483 if ( sopt->debug > 0 )
00484 fprintf(stderr, "-- (min_dist, max_dist) = (%f, %f)\n"
00485 " %d of %d nodes were out of bounds\n",
00486 min_dist, max_dist, oobc, N->ilen);
00487
00488 if ( final_computations(ddata, cdata, sopt, DSET_NVOX(p->gpar)) )
00489 RETURN(-1);
00490
00491 RETURN(0);
00492 }
00493
00494
00495
00496
00497
00498
00499 int final_computations(double *ddata, int *cdata, s2v_opts_t *sopt, int nvox)
00500 {
00501 int index;
00502
00503 ENTRY("final_computations");
00504
00505 if ( (sopt->debug > 1) && (sopt->dvox >= 0) && (sopt->dvox < nvox) )
00506 fprintf(stderr,"++ final: voxel %d, from values (%d, %f) ",
00507 sopt->dvox,cdata[sopt->dvox],ddata[sopt->dvox]);
00508
00509 switch ( sopt->map )
00510 {
00511 default:
00512 fprintf(stderr,"** fc: mapping %d not ready\n", sopt->map );
00513 RETURN(-1);
00514
00515 case E_SMAP_AVE:
00516
00517 for ( index = 0; index < nvox; index++ )
00518 if ( cdata[index] > 0 )
00519 ddata[index] /= cdata[index];
00520 break;
00521
00522 case E_SMAP_COUNT:
00523 case E_SMAP_MAX_ABS:
00524 case E_SMAP_MAX:
00525 case E_SMAP_MIN:
00526 case E_SMAP_MASK:
00527 case E_SMAP_MASK2:
00528 break;
00529 }
00530
00531 if ( (sopt->debug > 1) && (sopt->dvox >= 0) && (sopt->dvox < nvox) )
00532 fprintf(stderr,"to values (%d, %f)\n",
00533 cdata[sopt->dvox],ddata[sopt->dvox]);
00534
00535 RETURN(0);
00536 }
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552 int insert_list( node_list_t * N, param_t * p, s2v_opts_t * sopt,
00553 THD_fvec3 * pary, int nindex, double * ddata, int * cdata )
00554 {
00555 THD_ivec3 i3;
00556 int vindex, prev_vind;
00557 int step, node;
00558 int nx, ny, debug;
00559
00560 ENTRY("insert_list");
00561
00562 if ( !N || !p || !sopt || !pary || nindex < 0 || !ddata || !cdata )
00563 {
00564 fprintf(stderr,"** ID params (%p,%p,%p,%p,%d,%p,%p)\n",
00565 N, p, sopt, pary, nindex, ddata, cdata);
00566 RETURN(-1);
00567 }
00568
00569 node = N->ilist[nindex];
00570
00571 debug = sopt->debug && node == sopt->dnode;
00572
00573 if ( debug )
00574 fprintf(stderr,"++ point list for N[%d] = %d:\n", nindex, node);
00575
00576 nx = DSET_NX(p->gpar);
00577 ny = DSET_NY(p->gpar);
00578 prev_vind = -1;
00579 for ( step = 0; step < sopt->f_steps; step++ )
00580 {
00581 if ( f3mm_out_of_bounds( pary+step, &p->f3mm_min, &p->f3mm_max ) )
00582 continue;
00583
00584 i3 = THD_3dmm_to_3dind(p->gpar, pary[step]);
00585 vindex = i3.ijk[0] + nx * (i3.ijk[1] + ny * i3.ijk[2]);
00586
00587 if ( debug )
00588 fprintf(stderr," %3d, vox %d, [%3d,%3d,%3d]: %f, %f, %f",
00589 step, vindex, i3.ijk[0], i3.ijk[1], i3.ijk[2],
00590 pary[step].xyz[0], pary[step].xyz[1], pary[step].xyz[2]);
00591
00592 if ( sopt->cmask && !sopt->cmask[vindex] )
00593 {
00594 if ( debug ) fprintf(stderr, " : skip (mask)\n");
00595 continue;
00596 }
00597
00598 if ( (vindex == prev_vind) && (sopt->f_index == S2V_F_INDEX_VOXEL) )
00599 {
00600 if ( debug ) fprintf(stderr, " : skip (repeat voxel)\n");
00601 continue;
00602 }
00603
00604
00605
00606 prev_vind = vindex;
00607
00608 if ( debug )
00609 fprintf( stderr, " : (old) %d %f", cdata[vindex],ddata[vindex]);
00610
00611 if (insert_value(sopt, ddata, cdata, vindex, node, N->fdata[nindex]))
00612 RETURN(-1);
00613
00614 if ( debug )
00615 fprintf( stderr, " : (new) %d %f\n", cdata[vindex],ddata[vindex]);
00616 }
00617
00618 RETURN(0);
00619 }
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629 int insert_value(s2v_opts_t * sopt, double *dv, int *cv, int vox, int node,
00630 float value)
00631 {
00632 ENTRY("insert_value");
00633
00634 if ( !sopt || !dv || !cv || vox < 0 )
00635 {
00636 fprintf(stderr, "** IV: bad params (%p,%p,%p,%d)\n",sopt,dv,cv,vox);
00637 RETURN(-1);
00638 }
00639
00640 if ( (sopt->debug > 1) && (vox == sopt->dvox) )
00641 fprintf(stderr,"++ voxel %d, node %d, values (%d, %f) ",
00642 vox, node, cv[vox], dv[vox]);
00643
00644 switch ( sopt->map )
00645 {
00646 default:
00647 fprintf(stderr,"** IV: mapping %d not ready\n", sopt->map );
00648 RETURN(-1);
00649
00650 case E_SMAP_AVE:
00651 if ( cv[vox] == 0 )
00652 dv[vox] = value;
00653 else
00654 dv[vox] += value;
00655 break;
00656
00657 case E_SMAP_COUNT:
00658 if ( cv[vox] == 0 )
00659 dv[vox] = 1.0;
00660 else
00661 dv[vox]++;
00662 break;
00663
00664 case E_SMAP_MAX:
00665 if ( cv[vox] == 0 )
00666 dv[vox] = value;
00667 else if (value > dv[vox])
00668 dv[vox] = value;
00669
00670 break;
00671
00672 case E_SMAP_MAX_ABS:
00673 if ( cv[vox] == 0 )
00674 dv[vox] = value;
00675 else if (fabs(value) > fabs(dv[vox]))
00676 dv[vox] = value;
00677
00678 break;
00679
00680 case E_SMAP_MIN:
00681 if ( cv[vox] == 0 )
00682 dv[vox] = value;
00683 else if (value < dv[vox])
00684 dv[vox] = value;
00685
00686 break;
00687
00688 case E_SMAP_MASK:
00689 case E_SMAP_MASK2:
00690 dv[vox] = 1.0;
00691 break;
00692 }
00693
00694 cv[vox]++;
00695
00696 if ( (sopt->debug > 1) && (vox == sopt->dvox) )
00697 fprintf(stderr,"to (%d, %f)\n",cv[vox],dv[vox]);
00698
00699 RETURN(0);
00700 }
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710 int make_point_list( THD_fvec3 * list, THD_fvec3 * p1, THD_fvec3 * pn,
00711 s2v_opts_t * sopt )
00712 {
00713 float rat1, ratn;
00714 int step;
00715
00716 ENTRY("make_point_list");
00717
00718 if ( !list || !p1 || !pn || !sopt )
00719 {
00720 fprintf(stderr,"** mpl: bad params (%p,%p,%p,%p)\n", list,p1,pn,sopt);
00721 RETURN(-1);
00722 }
00723
00724 list[0] = *p1;
00725
00726 for ( step = 1; step < sopt->f_steps; step++ )
00727 {
00728 ratn = step / (sopt->f_steps - 1.0);
00729 rat1 = 1.0 - ratn;
00730
00731 list[step].xyz[0] = rat1 * p1->xyz[0] + ratn * pn->xyz[0];
00732 list[step].xyz[1] = rat1 * p1->xyz[1] + ratn * pn->xyz[1];
00733 list[step].xyz[2] = rat1 * p1->xyz[2] + ratn * pn->xyz[2];
00734 }
00735
00736 RETURN(0);
00737 }
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747 int adjust_endpts( s2v_opts_t * sopt, THD_fvec3 * p1, THD_fvec3 * pn )
00748 {
00749 THD_fvec3 f3_diff;
00750 float dist, factor;
00751
00752 ENTRY("adjust_endpts");
00753
00754 if ( !sopt || !p1 || !pn )
00755 {
00756 fprintf(stderr,"** ae: invalid params (%p,%p,%p)\n", sopt, p1, pn);
00757 RETURN(-1);
00758 }
00759
00760
00761 f3_diff.xyz[0] = pn->xyz[0] - p1->xyz[0];
00762 f3_diff.xyz[1] = pn->xyz[1] - p1->xyz[1];
00763 f3_diff.xyz[2] = pn->xyz[2] - p1->xyz[2];
00764
00765 dist = dist_f3mm( p1, pn );
00766
00767 if ( (sopt->f_p1_fr != 0.0) || (sopt->f_p1_mm != 0.0) )
00768 {
00769 if ( sopt->f_p1_fr != 0.0 )
00770 factor = sopt->f_p1_fr;
00771 else
00772 factor = (dist == 0.0) ? 0.0 : sopt->f_p1_mm / dist;
00773
00774 p1->xyz[0] += factor * f3_diff.xyz[0];
00775 p1->xyz[1] += factor * f3_diff.xyz[1];
00776 p1->xyz[2] += factor * f3_diff.xyz[2];
00777 }
00778
00779 if ( (sopt->f_pn_fr != 0.0) || (sopt->f_pn_mm != 0.0) )
00780 {
00781 if ( sopt->f_pn_fr != 0.0 )
00782 factor = sopt->f_pn_fr;
00783 else
00784 factor = (dist == 0.0) ? 0.0 : sopt->f_pn_mm / dist;
00785
00786 pn->xyz[0] += factor * f3_diff.xyz[0];
00787 pn->xyz[1] += factor * f3_diff.xyz[1];
00788 pn->xyz[2] += factor * f3_diff.xyz[2];
00789 }
00790
00791 switch ( sopt->map )
00792 {
00793 default:
00794 fprintf(stderr,"** ae: mapping %d not ready\n", sopt->map );
00795 RETURN(-1);
00796
00797 case E_SMAP_AVE:
00798 case E_SMAP_COUNT:
00799 case E_SMAP_MAX_ABS:
00800 case E_SMAP_MAX:
00801 case E_SMAP_MIN:
00802 case E_SMAP_MASK:
00803 case E_SMAP_MASK2:
00804 break;
00805 }
00806
00807 RETURN(0);
00808 }
00809
00810
00811
00812
00813
00814
00815
00816
00817 int set_node_list_data ( node_list_t *N, param_t *p, s2v_opts_t *sopt, int col )
00818 {
00819 float * fp, fval;
00820 int c, lposn;
00821
00822 ENTRY("set_node_list_data");
00823
00824 if ( !N || !p || !sopt || col < 0 )
00825 {
00826 fprintf(stderr,"** snld: bad params (%p,%p,%p,%d)\n", N, p, sopt, col);
00827 RETURN(-1);
00828 }
00829
00830 if ( sopt->debug > 1 )
00831 fprintf(stderr, "-- setting fdata for column %d\n", col);
00832
00833
00834 if ( !p->sdata_im )
00835 {
00836 fval = 1.0;
00837
00838
00839 if ( p->parser.pcode )
00840 {
00841 for ( c = 0; c < 26; c++ )
00842 p->parser.atoz[c] = 0.0;
00843 fval = PARSER_evaluate_one( p->parser.pcode, p->parser.atoz );
00844 }
00845
00846 for ( c = 0; c < N->ilen; c++ )
00847 N->fdata[c] = fval;
00848
00849 RETURN(0);
00850 }
00851
00852
00853 if ( col > (p->sdata_im->nx - 2) && !p->parser.pcode )
00854 {
00855 fprintf(stderr,"** snld error: col > nx-2 (%d > %d)\n",
00856 col, p->sdata_im->nx-2);
00857 RETURN(-1);
00858 }
00859 else if ( p->sdata_im->ny < N->ilen )
00860 {
00861 fprintf(stderr,"** snld error: ny < ilen (%d < %d)\n",
00862 p->sdata_im->ny, N->ilen);
00863 RETURN(-1);
00864 }
00865 else if ( !N->fdata )
00866 {
00867 fprintf(stderr,"** snld error: missing idata\n");
00868 RETURN(-1);
00869 }
00870 else if ( p->parser.pcode && col != 0 )
00871 {
00872 fprintf(stderr,"** snld error: cannot use parser with col = %d\n", col);
00873 RETURN(-1);
00874 }
00875
00876
00877
00878 fp = MRI_FLOAT_PTR( p->sdata_im ) + col+1;
00879
00880 for ( c = 0; c < N->ilen; c++ )
00881 {
00882 if ( p->parser.pcode )
00883 {
00884
00885 for ( lposn = 0; lposn < p->parser.max_sym; lposn++ )
00886 p->parser.atoz[lposn] = fp[lposn];
00887
00888 N->fdata[c] = PARSER_evaluate_one(p->parser.pcode, p->parser.atoz);
00889 }
00890 else
00891 N->fdata[c] = *fp;
00892
00893 fp += p->sdata_im->nx;
00894 }
00895
00896 RETURN(0);
00897 }
00898
00899
00900
00901
00902
00903
00904
00905
00906 int init_node_list (opts_t *opts, param_t *p, s2v_opts_t *sopt, node_list_t *N)
00907 {
00908 int nsurf, rv;
00909
00910 ENTRY("init_node_list");
00911
00912 if ( opts == NULL || p == NULL || sopt == NULL || N == NULL )
00913 {
00914 fprintf(stderr, "** inl - bad params (%p,%p,%p,%p)\n",opts,p,sopt,N);
00915 RETURN(-1);
00916 }
00917
00918 memset(N,0,sizeof(*N));
00919 N->nodes = NULL; N->fdata = NULL; N->ilist = NULL;
00920
00921 nsurf = SUMAg_N_DOv;
00922 if ( !p->sxyz_im && ((nsurf < 1) || (nsurf > S2V_MAX_SURFS)) )
00923 {
00924 fprintf(stderr,"** inl: SUMAg_N_DOv has invalid value of %d\n", nsurf);
00925 RETURN(-1);
00926 }
00927
00928 #if 0
00929
00930 if ( sopt->map == E_SMAP_MASK )
00931 nsurf = 1;
00932 else
00933 {
00934
00935 if ( ! p->sxyz_im && (SUMAg_N_DOv < 2) )
00936 {
00937 fprintf(stderr,"** function '%s' requires 2 surfaces\n",
00938 gs2v_map_names[sopt->map]);
00939 RETURN(-1);
00940 }
00941
00942 nsurf = 2;
00943 }
00944 #endif
00945
00946 if ( p->sxyz_im )
00947 rv = sxyz_1D_to_nlist( opts, sopt, p, N, &nsurf );
00948 else
00949 rv = surf_to_node_list( sopt, N, nsurf );
00950
00951 if ( sopt->debug > 1 && rv == 0 )
00952 fprintf(stderr,"++ node list initialized\n");
00953
00954 RETURN(rv);
00955 }
00956
00957
00958
00959
00960
00961
00962
00963 int sxyz_1D_to_nlist(opts_t * opts, s2v_opts_t * sopt, param_t * p,
00964 node_list_t * N, int * nsurf)
00965 {
00966 THD_fvec3 * fvp;
00967 float * fp;
00968 int c, sc;
00969
00970 ENTRY("sxyz_1D_to_nlist");
00971
00972 if ( !sopt || !p || !N )
00973 {
00974 fprintf(stderr,"** sxyz2nl: bad params (%p,%p,%p)\n",sopt,p,N);
00975 RETURN(-1);
00976 }
00977
00978 if ( !p->sxyz_im )
00979 {
00980 fprintf(stderr,"** missing sxyz_im for sxyz surf '%s'\n",
00981 opts->surf_xyz_1D_file);
00982 RETURN(-1);
00983 }
00984
00985 *nsurf = p->sxyz_im->nx / 3;
00986
00987 if ( p->sxyz_im->nx != 3 * *nsurf )
00988 {
00989 fprintf(stderr,"** sxyz surf '%s' has %d columns (%d expected)\n",
00990 opts->surf_xyz_1D_file, p->sxyz_im->nx, 3**nsurf);
00991 RETURN(-1);
00992 }
00993 else if ( p->sxyz_im->ny <= 0 )
00994 {
00995 fprintf(stderr,"** sxyz surf '%s': bad sxyz dimensions (%d,%d)\n",
00996 opts->surf_xyz_1D_file, p->sxyz_im->nx, p->sxyz_im->ny);
00997 RETURN(-1);
00998 }
00999
01000 N->depth = *nsurf;
01001 N->nnodes = p->sxyz_im->ny;
01002
01003 N->nodes = (THD_fvec3 *)malloc(N->depth*N->nnodes*sizeof(THD_fvec3));
01004 if ( N->nodes == NULL )
01005 {
01006 fprintf(stderr,"** failed to allocate %dx%d THD_fvec3's for nodes\n",
01007 N->nnodes, N->depth);
01008 RETURN(-1);
01009 }
01010
01011 fvp = N->nodes;
01012 for ( sc = 0; sc < *nsurf; sc++ )
01013 {
01014 fp = MRI_FLOAT_PTR( p->sxyz_im ) + sc * 3;
01015 for ( c = 0; c < N->nnodes; c++ )
01016 {
01017 fvp->xyz[0] = fp[0];
01018 fvp->xyz[1] = fp[1];
01019 fvp->xyz[2] = fp[2];
01020
01021 fp += p->sxyz_im->nx;
01022 fvp++;
01023 }
01024 }
01025
01026 if ( sopt->debug > 0 )
01027 {
01028 fprintf(stderr, "++ sxyz_1D nodes from '%s': nxyz = %d, nsurf = %d\n",
01029 opts->surf_xyz_1D_file, N->nnodes, N->depth);
01030 if ( sopt->dnode >= 0 && sopt->dnode <= p->sxyz_im->ny )
01031 {
01032 fprintf(stderr," debug node (%d) loc", sopt->dnode);
01033 for ( sc = 0; sc < *nsurf; sc++ )
01034 fprintf(stderr," : (%f, %f, %f)",
01035 N->nodes[sopt->dnode + sc*N->nnodes].xyz[0],
01036 N->nodes[sopt->dnode + sc*N->nnodes].xyz[1],
01037 N->nodes[sopt->dnode + sc*N->nnodes].xyz[2]);
01038 fputc('\n', stderr);
01039 }
01040 }
01041
01042 RETURN(0);
01043 }
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055 int surf_to_node_list ( s2v_opts_t * sopt, node_list_t * N, int nsurf )
01056 {
01057 SUMA_SurfaceObject ** so;
01058 THD_fvec3 * fvp;
01059 float * fp;
01060 int rv, nindex, sindex;
01061
01062 ENTRY("surf_to_node_list");
01063
01064 if ( sopt == NULL || N == NULL || nsurf < 0 || nsurf > 2 )
01065 {
01066 fprintf( stderr, "** anl: bad params (%p,%p,%d)\n",
01067 sopt, N, nsurf );
01068 RETURN(-1);
01069 }
01070
01071
01072 so = (SUMA_SurfaceObject **)calloc(nsurf, sizeof(SUMA_SurfaceObject *));
01073 if ( so == NULL )
01074 {
01075 fprintf( stderr, "** anl: failed to alloc %d surf pointers\n", nsurf );
01076 RETURN(-1);
01077 }
01078
01079 if ( (rv = get_mappable_surfs( so, nsurf, sopt->debug )) != nsurf )
01080 {
01081 fprintf( stderr, "** error: found %d (of %d) mappable surfaces\n",
01082 rv, nsurf );
01083 RETURN(-1);
01084 }
01085
01086
01087 N->depth = nsurf;
01088 N->nnodes = so[0]->N_Node;
01089 N->nodes = (THD_fvec3 *)malloc(N->depth * N->nnodes * sizeof(THD_fvec3));
01090 if ( N->nodes == NULL )
01091 {
01092 fprintf( stderr, "** cnlm: failed to allocate %d THD_fvec3 structs\n",
01093 N->depth * N->nnodes );
01094 free(so);
01095 RETURN(-1);
01096 }
01097
01098
01099
01100 fvp = N->nodes;
01101 for ( sindex = 0; sindex < N->depth; sindex++ )
01102 {
01103 if ( so[sindex]->N_Node != N->nnodes )
01104 {
01105 fprintf(stderr, "** surf #%d (%s) has %d nodes (but expected %d)\n",
01106 sindex,
01107 so[sindex]->Label ? so[sindex]->Label : "<unnamed>",
01108 so[sindex]->N_Node, N->nnodes );
01109 free( N->nodes ); N->nodes = NULL;
01110 free(so);
01111 RETURN(-1);
01112 }
01113
01114 for ( nindex = 0, fp = so[sindex]->NodeList;
01115 nindex < N->nnodes;
01116 nindex++, fp += 3 )
01117 {
01118 memcpy( fvp->xyz, fp, 3*sizeof(float) );
01119 fvp++;
01120 }
01121 }
01122
01123 if ( sopt->debug > 1 )
01124 fprintf( stderr, "++ allocated %d x %d (x %d) node list\n",
01125 N->depth, N->nnodes, (int)sizeof(THD_fvec3) );
01126
01127 free(so);
01128 RETURN(0);
01129 }
01130
01131
01132
01133
01134
01135
01136
01137
01138 int get_mappable_surfs( SUMA_SurfaceObject ** slist, int how_many, int debug )
01139 {
01140 SUMA_SurfaceObject * so;
01141 int count, socount = 0;
01142
01143 ENTRY("get_mappable_surfts");
01144
01145 if ( slist == NULL )
01146 {
01147 fprintf( stderr, "** gms: missing slist!\n" );
01148 RETURN(-1);
01149 }
01150
01151 for ( count = 0; count < SUMAg_N_DOv; count++ )
01152 {
01153 if ( ! SUMA_isSO(SUMAg_DOv[count]) )
01154 continue;
01155
01156 so = (SUMA_SurfaceObject *)SUMAg_DOv[count].OP;
01157
01158 if ( ! so->AnatCorrect )
01159 {
01160 if ( debug )
01161 fprintf(stderr,"-- surf #%d '%s', anat not correct, skipping\n",
01162 socount, CHECK_NULL_STR(so->Label));
01163 if ( debug > 1 )
01164 fprintf(stderr,"** consider adding the following to the "
01165 "surface definition in the spec file:\n"
01166 " Anatomical = Y\n");
01167 continue;
01168 }
01169
01170 if ( debug > 1 )
01171 {
01172 fprintf( stderr, "\n---------- surface #%d '%s' -----------\n",
01173 socount, CHECK_NULL_STR(so->Label) );
01174 SUMA_Print_Surface_Object( so, stderr );
01175 }
01176
01177 if ( socount < how_many )
01178 slist[socount] = so;
01179
01180 socount++;
01181 }
01182
01183 if ( debug > 1 )
01184 fprintf( stderr, "++ found %d mappable surfaces\n", socount );
01185
01186 RETURN(socount);
01187 }
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197 int set_smap_opts( opts_t * opts, param_t * p, s2v_opts_t * sopt )
01198 {
01199 int nsurf;
01200
01201 ENTRY("set_smap_opts");
01202
01203 memset( sopt, 0, sizeof(*sopt) );
01204 sopt->cmask = NULL;
01205
01206 if ( (sopt->map = check_map_func( opts->map_str )) == E_SMAP_INVALID )
01207 RETURN(-1);
01208
01209 sopt->datum = check_datum_type(opts->datum_str, DSET_BRICK_TYPE(p->gpar,0));
01210 if (sopt->datum < 0)
01211 RETURN(-1);
01212
01213 if ( opts->noscale == 1 )
01214 sopt->noscale = 1;
01215
01216 sopt->debug = opts->debug;
01217 sopt->dnode = opts->dnode;
01218 sopt->dvox = opts->dvox;
01219 sopt->sxyz_ori_gpar = opts->sxyz_ori_gpar;
01220 sopt->cmask = p->cmask;
01221
01222 if ( opts->f_steps < S2V_F_STEPS_MIN )
01223 {
01224
01225 for ( nsurf = 0; nsurf < S2V_MAX_SURFS && opts->snames[nsurf]; nsurf++ )
01226 ;
01227 if ( nsurf <= 0 )
01228 {
01229 fprintf(stderr,"** error: sso: no input surfaces\n");
01230 RETURN(-1);
01231 }
01232 sopt->f_steps = nsurf;
01233 }
01234 else
01235 sopt->f_steps = opts->f_steps;
01236
01237 sopt->f_index = S2V_F_INDEX_VOXEL;
01238 if ( (opts->f_index_str != NULL) &&
01239 ( !strncmp(opts->f_index_str, "point", 5) ||
01240 !strncmp(opts->f_index_str, "node", 4) ) )
01241 sopt->f_index = S2V_F_INDEX_POINT;
01242
01243 sopt->f_p1_fr = opts->f_p1_fr;
01244 sopt->f_pn_fr = opts->f_pn_fr;
01245 sopt->f_p1_mm = opts->f_p1_mm;
01246 sopt->f_pn_mm = opts->f_pn_mm;
01247
01248
01249 switch (sopt->map)
01250 {
01251 default:
01252
01253 case E_SMAP_AVE:
01254 case E_SMAP_MAX:
01255 case E_SMAP_MIN:
01256 case E_SMAP_MAX_ABS:
01257 break;
01258
01259 case E_SMAP_COUNT:
01260 case E_SMAP_MASK:
01261 case E_SMAP_MASK2:
01262 sopt->noscale = 1;
01263 break;
01264 }
01265
01266 if ( opts->debug > 1 )
01267 disp_s2v_opts_t( "++ s2v_opts_set :", sopt );
01268
01269 RETURN(0);
01270 }
01271
01272
01273
01274
01275
01276
01277 int final_clean_up ( node_list_t * N )
01278 {
01279 ENTRY("final_clean_up");
01280
01281 if ( ( SUMAg_DOv != NULL ) &&
01282 ( SUMA_Free_Displayable_Object_Vect(SUMAg_DOv, SUMAg_N_DOv) == 0 ) )
01283 fprintf(stderr, "** failed SUMA_Free_Displayable_Object_Vect()\n" );
01284
01285 if ( ( SUMAg_SVv != NULL ) &&
01286 ( SUMA_Free_SurfaceViewer_Struct_Vect(SUMAg_SVv, SUMAg_N_SVv) == 0 ) )
01287 fprintf( stderr, "** failed SUMA_Free_SurfaceViewer_Struct_Vect()\n" );
01288
01289 if ( ( SUMAg_CF != NULL ) && ( SUMA_Free_CommonFields(SUMAg_CF) == 0 ) )
01290 fprintf( stderr, "** failed SUMA_Free_CommonFields()\n" );
01291
01292 if ( N )
01293 {
01294 if ( N->nodes ) free( N->nodes );
01295 if ( N->fdata ) free( N->fdata );
01296 if ( N->ilist ) free( N->ilist );
01297 }
01298
01299 RETURN(0);
01300 }
01301
01302
01303
01304
01305
01306
01307 int read_surf_files ( opts_t * opts, param_t * p, SUMA_SurfSpecFile * spec,
01308 s2v_opts_t * sopt, node_list_t * N )
01309 {
01310 int rv;
01311
01312 ENTRY("read_surf_files");
01313
01314
01315 if ( opts->spec_file )
01316 {
01317 if ( (rv = fill_SUMA_structs(opts, spec)) != 0 )
01318 RETURN(rv);
01319 }
01320 else if ( opts->surf_xyz_1D_file )
01321 {
01322 if ( (rv = read_sxyz_1D( opts, p )) != 0 )
01323 RETURN(rv);
01324 }
01325 else
01326 {
01327 fprintf(stderr,"** missing spec or sdata file, exiting...\n");
01328 RETURN(-1);
01329 }
01330
01331
01332 if ( (rv = init_node_list(opts, p, sopt, N)) != 0 )
01333 RETURN(rv);
01334
01335
01336 if ( (rv = fill_node_list(opts, p, N)) != 0 )
01337 RETURN(rv);
01338
01339 if ( (rv = verify_parser_expr(opts, p)) != 0 )
01340 RETURN(rv);
01341
01342 RETURN(0);
01343 }
01344
01345
01346
01347
01348
01349
01350 int read_sxyz_1D ( opts_t * opts, param_t * p )
01351 {
01352 MRI_IMAGE * im;
01353
01354 ENTRY("read_sxyz_1D");
01355
01356 if ( !opts || !p )
01357 {
01358 fprintf(stderr,"** rsxyz1D: bad params (%p,%p)\n", opts, p);
01359 RETURN(-1);
01360 }
01361
01362 if ( (im = mri_read_1D(opts->surf_xyz_1D_file)) == NULL )
01363 {
01364 fprintf(stderr,"** failed to read sxyz 1D file '%s'\n",
01365 opts->surf_xyz_1D_file);
01366 RETURN(-1);
01367 }
01368 if ( im->nx < 1 || im->ny < 3 )
01369 {
01370 fprintf(stderr,"** bad sxyz file '%s'?\n", opts->surf_xyz_1D_file);
01371 RETURN(-1);
01372 }
01373
01374
01375 p->sxyz_im = mri_transpose(im);
01376 mri_free(im);
01377
01378 if ( opts->debug > 0 )
01379 fprintf(stderr,"++ read 1D xyz surface file '%s' (nx = %d, ny = %d)\n",
01380 opts->surf_xyz_1D_file, p->sxyz_im->nx, p->sxyz_im->ny );
01381
01382 RETURN(0);
01383 }
01384
01385
01386
01387
01388
01389
01390
01391
01392
01393
01394 int fill_node_list ( opts_t * opts, param_t * p, node_list_t * N )
01395 {
01396 int rv;
01397
01398 ENTRY("fill_node_list");
01399
01400 if ( !opts || !p || !N )
01401 {
01402 fprintf(stderr,"** fill_node_list: bad params (%p,%p,%p)\n",opts,p,N);
01403 RETURN(-1);
01404 }
01405
01406 p->nsubs = 1;
01407
01408 if ( opts->sdata_file_1D )
01409 {
01410 if ( (rv = sdata_from_1D( opts, p, N )) != 0 )
01411 RETURN(rv);
01412 if ( ! p->parser.pcode )
01413 p->nsubs = p->sdata_im->nx - 1;
01414 }
01415 else
01416 {
01417 if ( (rv = sdata_from_default( N )) != 0 )
01418 RETURN(rv);
01419 }
01420
01421 if ( (rv = verify_node_list( N )) != 0 )
01422 RETURN(rv);
01423
01424 if ( opts->debug > 1 )
01425 disp_node_list_t( "++ node list filled: ", N );
01426
01427 RETURN(0);
01428 }
01429
01430
01431
01432
01433
01434
01435 int verify_parser_expr ( opts_t * opts, param_t * p )
01436 {
01437 int max_used;
01438
01439 ENTRY("verify_parser_expr");
01440
01441 if ( !opts || !p )
01442 {
01443 fprintf(stderr,"** vpe: invalid params (%p,%p)\n", opts, p);
01444 RETURN(-1);
01445 }
01446
01447
01448 if ( ! p->parser.pcode )
01449 RETURN(0);
01450
01451 for ( max_used = 25; max_used >= 0; max_used-- )
01452 if ( p->parser.has_sym[max_used] )
01453 break;
01454 max_used++;
01455 p->parser.max_sym = max_used;
01456
01457
01458 if ( max_used > 0 )
01459 {
01460 if ( !p->sdata_im )
01461 {
01462 fprintf(stderr, "** parser expression requires surface data\n"
01463 " (see '-sdata_1D')\n");
01464 RETURN(-1);
01465 }
01466 else if ( max_used > p->sdata_im->nx - 1 )
01467 {
01468 fprintf(stderr,
01469 "** error: not enough surface values for expression\n"
01470 " svals = %d, exp_vals = %d, expr = '%s'\n",
01471 p->sdata_im->nx - 1, max_used, opts->data_expr);
01472 RETURN(-1);
01473 }
01474 }
01475
01476 if ( opts->debug > 1 )
01477 fprintf(stderr,"-- surf_vals = %d, expr_vals = %d\n",
01478 p->sdata_im ? (p->sdata_im->nx - 1) : 0, max_used);
01479
01480 RETURN(0);
01481 }
01482
01483
01484
01485
01486
01487
01488 int verify_node_list ( node_list_t * N )
01489 {
01490 int icount, errs = 0;
01491
01492 ENTRY("verify_node_list");
01493
01494 if ( !N )
01495 {
01496 fprintf(stderr, "** vnl - no node list\n" );
01497 RETURN(-1);
01498 }
01499
01500 if ( !N->nodes || !N->ilist )
01501 {
01502 fprintf(stderr,"** missing nodes or ilist\n" );
01503 errs++;
01504 }
01505
01506 if ( N->depth < 1 || N->nnodes < 1 || N->ilen < 1 )
01507 {
01508 fprintf(stderr,"** invalid depth, nnodes or ilen" );
01509 errs++;
01510 }
01511
01512 if ( errs )
01513 {
01514 disp_node_list_t("** invalid data : ", N );
01515 RETURN(-errs);
01516 }
01517
01518
01519 for ( icount = 0; icount < N->ilen; icount++ )
01520 if ( N->ilist[icount] < 0 || N->ilist[icount] >= N->nnodes )
01521 {
01522 fprintf(stderr,"** surf data index number %d is out of range:\n"
01523 " index = %d, range is [%d,%d]\n",
01524 icount, N->ilist[icount], 0, N->nnodes-1);
01525 RETURN(-10);
01526 }
01527
01528
01529 if ( (N->fdata = (float *)malloc(N->ilen*sizeof(float))) == NULL )
01530 {
01531 fprintf(stderr,"** vnl: failed to allocate %d floats\n",N->ilen);
01532 free(N->ilist);
01533 RETURN(-20);
01534 }
01535
01536 RETURN(0);
01537 }
01538
01539
01540
01541
01542
01543
01544 int sdata_from_default ( node_list_t * N )
01545 {
01546 int c;
01547
01548 ENTRY("sdata_from_default");
01549
01550 if ( !N )
01551 RETURN(-1);
01552
01553 if ( (N->ilist = (int *)malloc(N->nnodes * sizeof(int))) == NULL )
01554 {
01555 fprintf(stderr,"** failed to allocate %d ints for ilist\n",N->nnodes);
01556 RETURN(-1);
01557 }
01558
01559 N->ilen = N->nnodes;
01560
01561
01562 for ( c = 0; c < N->ilen; c++ )
01563 N->ilist[c] = c;
01564
01565 RETURN(0);
01566 }
01567
01568
01569
01570
01571
01572
01573 int sdata_from_1D ( opts_t * opts, param_t * p, node_list_t * N )
01574 {
01575 MRI_IMAGE * im;
01576 float * fim;
01577 int c;
01578
01579 ENTRY("sdata_from_1D");
01580
01581 if ( !opts || !N || !p )
01582 RETURN(-1);
01583
01584 if ( (im = mri_read_1D(opts->sdata_file_1D)) == NULL )
01585 {
01586 fprintf(stderr,"** failed to read 1D file '%s'\n", opts->sdata_file_1D);
01587 RETURN(-2);
01588 }
01589
01590 p->sdata_im = mri_transpose(im);
01591 mri_free(im);
01592
01593 if ( p->sdata_im->nx < 2 || p->sdata_im->ny < 1 )
01594 {
01595 fprintf(stderr,"** bad (%d x %d) surf data 1D file '%s'\n",
01596 p->sdata_im->ny, p->sdata_im->nx, opts->sdata_file_1D);
01597 RETURN(-3);
01598 }
01599
01600 N->ilen = p->sdata_im->ny;
01601
01602 if ( opts->debug > 0 )
01603 fprintf(stderr,"++ read 1D surface file '%s' (nx = %d, ny = %d)\n",
01604 opts->sdata_file_1D, p->sdata_im->nx, p->sdata_im->ny );
01605
01606
01607 if ( (N->ilist = (int *)malloc(N->ilen*sizeof(int))) == NULL )
01608 {
01609 fprintf(stderr,"** sf1D a: failed to allocate %d ints\n", N->ilen);
01610 RETURN(-1);
01611 }
01612
01613
01614 fim = MRI_FLOAT_PTR( p->sdata_im );
01615 for ( c = 0; c < N->ilen; c++, fim += p->sdata_im->nx )
01616 N->ilist[c] = (int)*fim;
01617
01618 RETURN(0);
01619 }
01620
01621
01622
01623
01624
01625
01626
01627
01628 int fill_SUMA_structs ( opts_t * opts, SUMA_SurfSpecFile * spec )
01629 {
01630 int debug = 0, rv;
01631 ENTRY("fill_SUMA_structs");
01632
01633 if ( opts->debug > 2 )
01634 debug = 1;
01635
01636 if ( debug )
01637 fputs( "-- SUMA_Create_CommonFields()...\n", stderr );
01638
01639 if ( opts->spec_file == NULL )
01640 RETURN(-1);
01641
01642
01643 SUMAg_CF = SUMA_Create_CommonFields();
01644
01645 if ( SUMAg_CF == NULL )
01646 {
01647 fprintf( stderr, "** failed SUMA_Create_CommonFields(), exiting...\n" );
01648 RETURN(-1);
01649 }
01650
01651
01652 if ( opts->debug > 3 )
01653 {
01654 SUMAg_CF->MemTrace = 1;
01655
01656 if ( opts->debug > 4 )
01657 SUMAg_CF->InOut_Notify = 1;
01658 }
01659
01660 if ( debug )
01661 fputs( "-- SUMA_Alloc_DisplayObject_Struct()...\n", stderr );
01662
01663 SUMAg_DOv = SUMA_Alloc_DisplayObject_Struct(SUMA_MAX_DISPLAYABLE_OBJECTS);
01664
01665 if ( debug )
01666 fputs( "-- SUMA_Read_SpecFile()...\n", stderr );
01667
01668 if ( SUMA_Read_SpecFile( opts->spec_file, spec) == 0 )
01669 {
01670 fprintf( stderr, "** failed SUMA_Read_SpecFile(), exiting...\n" );
01671 RETURN(-1);
01672 }
01673
01674 if ( debug )
01675 SUMA_ShowSpecStruct(spec, stderr, 1);
01676
01677 rv = SUMA_spec_select_surfs(spec, opts->snames, S2V_MAX_SURFS, opts->debug);
01678 if ( rv < 1 )
01679 {
01680 if ( rv == 0 )
01681 fprintf(stderr,"** no named surfaces found in spec file\n");
01682 RETURN(-1);
01683 }
01684
01685 if ( opts->debug > 0 )
01686 SUMA_ShowSpecStruct(spec, stderr, 1);
01687
01688 if ( SUMA_spec_set_map_refs(spec, opts->debug) != 0 )
01689 RETURN(-1);
01690
01691
01692 if ( spec->N_Groups != 1 )
01693 {
01694 fprintf( stderr,"** error: N_Groups <%d> must be 1 in spec file <%s>\n",
01695 spec->N_Groups, opts->spec_file );
01696 RETURN(-1);
01697 }
01698
01699 if ( debug )
01700 fputs( "-- SUMA_LoadSpec_eng()...\n", stderr );
01701
01702
01703 if (SUMA_LoadSpec_eng(spec,SUMAg_DOv,&SUMAg_N_DOv,opts->sv_file,debug,
01704 SUMAg_CF->DsetList) == 0)
01705 {
01706 fprintf( stderr, "** error: failed SUMA_LoadSpec_eng(), exiting...\n" );
01707 RETURN(-1);
01708 }
01709
01710 if ( opts->debug > 1 )
01711 fprintf(stderr, "++ %d surfaces loaded.\n", spec->N_Surfs );
01712
01713 RETURN(0);
01714 }
01715
01716
01717
01718
01719
01720
01721 int init_options ( opts_t * opts, int argc, char * argv [] )
01722 {
01723 int ac, ind;
01724
01725 ENTRY("init_options");
01726
01727 if ( argc < 2 )
01728 {
01729 usage( PROG_NAME, S2V_USE_LONG );
01730 RETURN(-1);
01731 }
01732
01733
01734 memset( opts, 0, sizeof( opts_t) );
01735 opts->gpar_file = NULL; opts->oset_file = NULL;
01736 opts->spec_file = NULL; opts->sv_file = NULL;
01737 opts->surf_xyz_1D_file = NULL; opts->sdata_file_1D = NULL;
01738 opts->sdata_file_niml = NULL; opts->cmask_cmd = NULL;
01739 opts->data_expr = NULL; opts->map_str = NULL;
01740 opts->datum_str = NULL; opts->f_index_str = NULL;
01741 opts->snames[0] = NULL; opts->snames[1] = NULL;
01742
01743
01744 opts->dnode = -1;
01745 opts->dvox = -1;
01746
01747 for ( ac = 1; ac < argc; ac++ )
01748 {
01749
01750 if ( ! strncmp(argv[ac], "-cmask", 6) )
01751 {
01752 if ( (ac+1) >= argc )
01753 {
01754 fputs( "option usage: -cmask COMMAND\n\n", stderr );
01755 usage( PROG_NAME, S2V_USE_SHORT );
01756 RETURN(-1);
01757 }
01758
01759 opts->cmask_cmd = argv[++ac];
01760 }
01761 else if ( ! strncmp(argv[ac], "-data_expr", 8) )
01762 {
01763 if ( (ac+1) >= argc )
01764 {
01765 fputs( "option usage: -data_expr EXPR\n\n", stderr );
01766 usage( PROG_NAME, S2V_USE_SHORT );
01767 RETURN(-1);
01768 }
01769
01770 opts->data_expr = argv[++ac];
01771 }
01772 else if ( ! strncmp(argv[ac], "-datum", 6) )
01773 {
01774 if ( (ac+1) >= argc )
01775 {
01776 fputs( "option usage: -datum DTYPE\n\n", stderr );
01777 usage( PROG_NAME, S2V_USE_SHORT );
01778 RETURN(-1);
01779 }
01780
01781 opts->datum_str = argv[++ac];
01782 }
01783 else if ( ! strncmp(argv[ac], "-debug", 6) )
01784 {
01785 if ( (ac+1) >= argc )
01786 {
01787 fputs( "option usage: -debug LEVEL\n\n", stderr );
01788 usage( PROG_NAME, S2V_USE_SHORT );
01789 RETURN(-1);
01790 }
01791
01792 opts->debug = atoi(argv[++ac]);
01793 if ( opts->debug < 0 || opts->debug > S2V_DEBUG_MAX_LEV )
01794 {
01795 fprintf( stderr, "bad debug level <%d>, should be in [0,%d]\n",
01796 opts->debug, S2V_DEBUG_MAX_LEV );
01797 usage( PROG_NAME, S2V_USE_SHORT );
01798 RETURN(-1);
01799 }
01800 }
01801 else if ( ! strncmp(argv[ac], "-dnode", 6) )
01802 {
01803 if ( (ac+1) >= argc )
01804 {
01805 fputs( "option usage: -dnode DEBUG_NODE\n\n", stderr );
01806 usage( PROG_NAME, S2V_USE_SHORT );
01807 RETURN(-1);
01808 }
01809
01810 opts->dnode = atoi(argv[++ac]);
01811 }
01812 else if ( ! strncmp(argv[ac], "-dvoxel", 5) )
01813 {
01814 if ( (ac+1) >= argc )
01815 {
01816 fputs( "option usage: -dvoxel DEBUG_VOXEL\n\n", stderr );
01817 usage( PROG_NAME, S2V_USE_SHORT );
01818 RETURN(-1);
01819 }
01820
01821 opts->dvox = atoi(argv[++ac]);
01822 }
01823 else if ( ! strncmp(argv[ac], "-f_index", 7) )
01824 {
01825 if ( (ac+1) >= argc )
01826 {
01827 fputs( "option usage: -f_index INDEX_TYPE\n\n", stderr );
01828 usage( PROG_NAME, S2V_USE_SHORT );
01829 RETURN(-1);
01830 }
01831 opts->f_index_str = argv[++ac];
01832 }
01833 else if ( ! strncmp(argv[ac], "-f_p1_fr", 9) )
01834 {
01835 if ( (ac+1) >= argc )
01836 {
01837 fputs( "option usage: -f_p1_fr FRACTION\n\n", stderr );
01838 usage( PROG_NAME, S2V_USE_SHORT );
01839 RETURN(-1);
01840 }
01841
01842 opts->f_p1_fr = atof(argv[++ac]);
01843 }
01844 else if ( ! strncmp(argv[ac], "-f_pn_fr", 9) )
01845 {
01846 if ( (ac+1) >= argc )
01847 {
01848 fputs( "option usage: -f_pn_fr FRACTION\n\n", stderr );
01849 usage( PROG_NAME, S2V_USE_SHORT );
01850 RETURN(-1);
01851 }
01852
01853 opts->f_pn_fr = atof(argv[++ac]);
01854 }
01855 else if ( ! strncmp(argv[ac], "-f_p1_mm", 9) )
01856 {
01857 if ( (ac+1) >= argc )
01858 {
01859 fputs( "option usage: -f_p1_mm DISTANCE\n\n", stderr );
01860 usage( PROG_NAME, S2V_USE_SHORT );
01861 RETURN(-1);
01862 }
01863
01864 opts->f_p1_mm = atof(argv[++ac]);
01865 }
01866 else if ( ! strncmp(argv[ac], "-f_pn_mm", 9) )
01867 {
01868 if ( (ac+1) >= argc )
01869 {
01870 fputs( "option usage: -f_pn_mm DISTANCE\n\n", stderr );
01871 usage( PROG_NAME, S2V_USE_SHORT );
01872 RETURN(-1);
01873 }
01874
01875 opts->f_pn_mm = atof(argv[++ac]);
01876 }
01877 else if ( ! strncmp(argv[ac], "-f_steps", 6) )
01878 {
01879 if ( (ac+1) >= argc )
01880 {
01881 fputs( "option usage: -f_steps NUM_STEPS\n\n", stderr );
01882 usage( PROG_NAME, S2V_USE_SHORT );
01883 RETURN(-1);
01884 }
01885
01886 opts->f_steps = atoi(argv[++ac]);
01887 }
01888 else if ( ! strncmp(argv[ac], "-grid_parent", 5)||
01889 ! strncmp(argv[ac], "-inset", 6) ||
01890 ! strncmp(argv[ac], "-input", 6) )
01891
01892 {
01893 if ( (ac+1) >= argc )
01894 {
01895 fputs( "option usage: -grid_parent INPUT_DSET\n\n", stderr );
01896 usage( PROG_NAME, S2V_USE_SHORT );
01897 RETURN(-1);
01898 }
01899
01900 opts->gpar_file = argv[++ac];
01901 }
01902 else if ( ! strncmp(argv[ac], "-help", 5) )
01903 {
01904 usage( PROG_NAME, S2V_USE_LONG );
01905 RETURN(-1);
01906 }
01907 else if ( ! strncmp(argv[ac], "-hist", 5) )
01908 {
01909 usage( PROG_NAME, S2V_USE_HIST );
01910 RETURN(-1);
01911 }
01912 else if ( ! strncmp(argv[ac], "-map_func", 4) )
01913 {
01914 if ( (ac+1) >= argc )
01915 {
01916 fputs( "option usage: -map_func FUNCTION\n\n", stderr );
01917 RETURN(-1);
01918 }
01919
01920 opts->map_str = argv[++ac];
01921 }
01922 else if ( ! strncmp(argv[ac], "-noscale", 4) )
01923 {
01924 opts->noscale = 1;
01925 }
01926 else if ( ! strncmp(argv[ac], "-prefix", 4) )
01927 {
01928 if ( (ac+1) >= argc )
01929 {
01930 fputs( "option usage: -prefix OUTPUT_PREFIX\n\n", stderr );
01931 usage( PROG_NAME, S2V_USE_SHORT );
01932 RETURN(-1);
01933 }
01934
01935 opts->oset_file = argv[++ac];
01936 }
01937 else if ( ! strncmp(argv[ac], "-sdata_1D", 9) )
01938 {
01939 if ( (ac+1) >= argc )
01940 {
01941 fputs( "option usage: -sdata_1D SURF_DATA.1D\n\n", stderr );
01942 usage( PROG_NAME, S2V_USE_SHORT );
01943 RETURN(-1);
01944 }
01945
01946 opts->sdata_file_1D = argv[++ac];
01947 }
01948 else if ( ! strncmp(argv[ac], "-sdata_niml", 9) )
01949 {
01950 if ( (ac+1) >= argc )
01951 {
01952 fputs( "option usage: -sdata_niml SURF_DATA.niml\n\n", stderr );
01953 usage( PROG_NAME, S2V_USE_SHORT );
01954 RETURN(-1);
01955 }
01956
01957 opts->sdata_file_niml = argv[++ac];
01958 }
01959 else if ( ! strncmp(argv[ac], "-spec", 5) )
01960 {
01961 if ( (ac+1) >= argc )
01962 {
01963 fputs( "option usage: -spec SPEC_FILE\n\n", stderr );
01964 usage( PROG_NAME, S2V_USE_SHORT );
01965 RETURN(-1);
01966 }
01967
01968 opts->spec_file = argv[++ac];
01969 }
01970 else if ( ! strncmp(argv[ac], "-surf_xyz_1D", 9) )
01971 {
01972 if ( (ac+1) >= argc )
01973 {
01974 fputs( "option usage: -surf_xyz_1D NODE_FILE\n\n", stderr );
01975 usage( PROG_NAME, S2V_USE_SHORT );
01976 RETURN(-1);
01977 }
01978
01979 opts->surf_xyz_1D_file = argv[++ac];
01980 }
01981 else if ( ! strncmp(argv[ac], "-surf_", 6) )
01982 {
01983 if ( (ac+1) >= argc )
01984 {
01985 fputs( "option usage: -surf_X SURF_NAME\n\n", stderr );
01986 usage( PROG_NAME, S2V_USE_SHORT );
01987 RETURN(-1);
01988 }
01989 ind = argv[ac][6] - 'A';
01990 if ( (ind < 0) || (ind >= S2V_MAX_SURFS) )
01991 {
01992 fprintf(stderr,"** -surf_X option: '%s' out of range,\n"
01993 " use one of '-surf_A' through '-surf_%c'\n",
01994 argv[ac], 'A'+S2V_MAX_SURFS-1);
01995 RETURN(-1);
01996 }
01997
01998 opts->snames[ind] = argv[++ac];
01999 }
02000 else if ( ! strncmp(argv[ac], "-sv", 3) )
02001 {
02002 if ( (ac+1) >= argc )
02003 {
02004 fputs( "option usage: -sv SURFACE_VOLUME\n\n", stderr );
02005 usage( PROG_NAME, S2V_USE_SHORT );
02006 RETURN(-1);
02007 }
02008
02009 opts->sv_file = argv[++ac];
02010 }
02011 else if ( ! strncmp(argv[ac], "-sxyz_orient_as_gpar", 20) )
02012 {
02013 opts->sxyz_ori_gpar = 1;
02014 }
02015 else if ( ! strncmp(argv[ac], "-version", 2) )
02016 {
02017 usage( PROG_NAME, S2V_USE_VERSION );
02018 RETURN(-1);
02019 }
02020 else
02021 {
02022 fprintf( stderr, "invalid option <%s>\n", argv[ac] );
02023 usage( PROG_NAME, S2V_USE_SHORT );
02024 RETURN(-1);
02025 }
02026 }
02027
02028 if ( opts->debug > 1 )
02029 disp_opts_t ( "++ opts read: ", opts );
02030
02031 RETURN(0);
02032 }
02033
02034
02035
02036
02037
02038
02039
02040
02041 int validate_options ( opts_t * opts, param_t * p )
02042 {
02043 ENTRY("validate_options");
02044
02045
02046 memset( p, 0, sizeof(*p) );
02047 p->gpar = NULL; p->oset = NULL;
02048 p->sxyz_im = NULL; p->sdata_im = NULL;
02049 p->parser.pcode = NULL; p->cmask = NULL;
02050
02051 if ( check_map_func( opts->map_str ) == E_SMAP_INVALID )
02052 RETURN(-1);
02053
02054 if ( validate_datasets( opts, p ) != 0 )
02055 RETURN(-1);
02056
02057 if ( validate_surface( opts, p ) != 0 )
02058 RETURN(-1);
02059
02060 if ( opts->data_expr )
02061 {
02062 PARSER_set_printout(1);
02063 p->parser.pcode = PARSER_generate_code( opts->data_expr );
02064 if ( p->parser.pcode == NULL )
02065 {
02066 fprintf(stderr,"** failed to generate code from expression '%s'\n",
02067 opts->data_expr);
02068 RETURN(-1);
02069 }
02070 PARSER_mark_symbols( p->parser.pcode, p->parser.has_sym );
02071
02072 if ( opts->debug > 1 )
02073 disp_parser_t( "-- PARSER expr okay: ", &p->parser );
02074 }
02075
02076 if ( opts->debug > 1 )
02077 disp_param_t( "++ params set: ", p );
02078
02079 RETURN(0);
02080 }
02081
02082
02083
02084
02085
02086
02087
02088
02089
02090
02091
02092 int check_datum_type ( char * datum_str, int default_type )
02093 {
02094 int c, dt = -1;
02095
02096 ENTRY("check_datum_type");
02097
02098 if ( datum_str )
02099 {
02100 for ( c = 0; c <= MRI_rgba; c++ )
02101 if ( ! strcmp( datum_str, MRI_TYPE_name[c] ) )
02102 {
02103 dt = c;
02104 break;
02105 }
02106
02107
02108 if ( c > MRI_rgba )
02109 {
02110 fprintf( stderr, "** invalid datum type name '%s'\n",
02111 datum_str );
02112 RETURN(-1);
02113 }
02114 }
02115 else
02116 {
02117 dt = default_type;
02118 }
02119
02120 if ( ( dt != MRI_byte ) &&
02121 ( dt != MRI_short ) &&
02122 ( dt != MRI_float ) &&
02123 ( dt != MRI_double ) )
02124 {
02125 fprintf( stderr, "** data type '%s' is not supported\n",
02126 MRI_TYPE_name[dt] );
02127 RETURN(-1);
02128 }
02129
02130 RETURN(dt);
02131 }
02132
02133
02134
02135
02136
02137
02138
02139
02140
02141 int check_map_func ( char * map_str )
02142 {
02143 int map;
02144
02145 ENTRY("check_map_func");
02146
02147 if ( map_str == NULL )
02148 {
02149 fprintf( stderr, "** missing option: '-map_func FUNCTION'\n" );
02150 RETURN(E_SMAP_INVALID);
02151 }
02152
02153 map = s2v_map_type( map_str );
02154
02155 if ( map == E_SMAP_INVALID )
02156 {
02157 fprintf( stderr, "** invalid map string '%s'\n", map_str );
02158 RETURN(-1);
02159 }
02160
02161 RETURN(map);
02162 }
02163
02164
02165
02166
02167
02168
02169
02170
02171
02172 int s2v_map_type ( char * map_str )
02173 {
02174 s2v_map_num map;
02175
02176 ENTRY("s2v_map_type");
02177
02178 if ( map_str == NULL )
02179 {
02180 fprintf( stderr, "** s2v_map_type: missing map_str parameter\n" );
02181 RETURN((int)E_SMAP_INVALID);
02182 }
02183
02184 if ( sizeof(gs2v_map_names) / sizeof(char *) != (int)E_SMAP_FINAL )
02185 {
02186 fprintf( stderr, "** error: gs2v_map_names/s2v_map_num mis-match\n");
02187 RETURN((int)E_SMAP_INVALID);
02188 }
02189
02190 for ( map = E_SMAP_MASK; map < E_SMAP_FINAL; map++ )
02191 if ( !strcmp( map_str, gs2v_map_names[map] ) )
02192 RETURN((int)map);
02193
02194 RETURN((int)E_SMAP_INVALID);
02195 }
02196
02197
02198
02199
02200
02201
02202
02203
02204
02205 int validate_surface ( opts_t * opts, param_t * p )
02206 {
02207 int errs = 0;
02208
02209 ENTRY("validate_surface");
02210
02211 if ( ! opts->surf_xyz_1D_file )
02212 {
02213 if ( opts->spec_file == NULL )
02214 {
02215 fprintf( stderr, "** missing '-spec_file SPEC_FILE' option\n" );
02216 errs++;
02217 }
02218
02219 if ( opts->sv_file == NULL )
02220 {
02221 fprintf( stderr, "** missing '-sv SURF_VOL' option\n" );
02222 errs++;
02223 }
02224
02225 if ( opts->snames[0] == NULL )
02226 {
02227 fprintf( stderr, "** missing '-surf_A SURF_NAME' option\n" );
02228 errs++;
02229 }
02230 }
02231 else if ( opts->spec_file != NULL )
02232 {
02233 fprintf(stderr,
02234 "** cannot use both spec and xyz surface files, %s and %s\n",
02235 opts->spec_file, opts->surf_xyz_1D_file);
02236 errs++;
02237 }
02238
02239 if ( opts->sdata_file_1D && opts->sdata_file_niml )
02240 {
02241 fprintf(stderr,"** cannot use both NIML and 1D surface files\n");
02242 errs++;
02243 }
02244
02245
02246 if ( opts->sdata_file_niml )
02247 {
02248 fprintf(stderr,"** sorry, the niml feature is coming soon...\n");
02249 errs++;
02250 }
02251
02252 if ( errs > 0 )
02253 RETURN(-1);
02254
02255 RETURN(0);
02256 }
02257
02258
02259
02260
02261
02262
02263
02264
02265
02266
02267
02268
02269 int validate_datasets( opts_t * opts, param_t * p )
02270 {
02271 ENTRY("validate_datasets");
02272
02273 p->gpar = THD_open_dataset( opts->gpar_file );
02274
02275 if ( !ISVALID_DSET(p->gpar) )
02276 {
02277 if ( opts->gpar_file == NULL )
02278 fprintf( stderr, "** error: missing '-grid_parent DSET' option\n" );
02279 else
02280 fprintf( stderr, "** error: invalid input dataset '%s'\n",
02281 opts->gpar_file);
02282 RETURN(-1);
02283 }
02284 else if ( DSET_BRICK_TYPE(p->gpar, 0) == MRI_complex )
02285 {
02286 fprintf(stderr,
02287 "** failure: cannot deal with complex-valued dataset, '%s'\n",
02288 opts->gpar_file);
02289 RETURN(-1);
02290 }
02291
02292 p->nvox = DSET_NVOX( p->gpar );
02293 set_3dmm_bounds( p->gpar, &p->f3mm_min, &p->f3mm_max );
02294
02295 if ( ! THD_filename_ok( opts->oset_file ) )
02296 {
02297 fprintf( stderr, "** illegal output prefix: '%s'\n",
02298 opts->oset_file ? opts->oset_file : "<none>" );
02299 RETURN(-1);
02300 }
02301
02302
02303
02304
02305 if ( opts->cmask_cmd != NULL )
02306 {
02307 int clen = strlen( opts->cmask_cmd );
02308 char * cmd;
02309
02310
02311 cmd = (char *)malloc((clen + 1) * sizeof(char));
02312 strcpy( cmd, opts->cmask_cmd );
02313
02314 p->cmask = EDT_calcmask( cmd, &p->ncmask );
02315
02316 free( cmd );
02317
02318 if ( p->cmask == NULL )
02319 {
02320 fprintf( stderr, "** failure: cannot compute mask from option:\n"
02321 " -cmask '%s'\n", opts->cmask_cmd );
02322 RETURN(-1);
02323 }
02324 if ( p->ncmask != p->nvox )
02325 {
02326 fprintf( stderr, "** error: input and cmask datasets do not have "
02327 "the same dimensions\n" );
02328 RETURN(-1);
02329 }
02330 if ( ( p->ccount = THD_countmask( p->ncmask, p->cmask ) ) <= 0 )
02331 {
02332 fprintf( stderr, "** Warning! No voxels in computed cmask!\n" );
02333
02334 }
02335 }
02336
02337 if ( opts->debug > 0 )
02338 {
02339 fprintf( stderr, "++ input dset has nvox = %d, nvals = %d",
02340 p->nvox, DSET_NVALS(p->gpar) );
02341 if ( p->cmask == NULL )
02342 fputc( '\n', stderr );
02343 else
02344 fprintf( stderr, " (%d voxels in mask)\n", p->ccount );
02345 }
02346
02347 RETURN(0);
02348 }
02349
02350
02351
02352
02353
02354
02355
02356
02357
02358 int usage ( char * prog, int level )
02359 {
02360 ENTRY("usage");
02361
02362 if ( level == S2V_USE_SHORT )
02363 fprintf(stderr,
02364 " usage: %s [options] -spec SPEC_FILE -surf_A SURF_NAME \\\n"
02365 " -grid_parent AFNI_DSET -sv SURF_VOL \\\n"
02366 " -map_func MAP_FUNC -prefix OUTPUT_DSET\n"
02367 "usage: %s -help\n",
02368 prog, prog );
02369 else if ( level == S2V_USE_LONG )
02370 {
02371 printf(
02372 "\n"
02373 "%s - map data from a surface domain to an AFNI volume domain\n"
02374 "\n"
02375 " usage: %s [options] -spec SPEC_FILE -surf_A SURF_NAME \\\n"
02376 " -grid_parent AFNI_DSET -sv SURF_VOL \\\n"
02377 " -map_func MAP_FUNC -prefix OUTPUT_DSET\n"
02378 "\n"
02379 " This program is meant to take as input a pair of surfaces,\n"
02380 " optionally including surface data, and an AFNI grid parent\n"
02381 " dataset, and to output a new AFNI dataset consisting of the\n"
02382 " surface data mapped to the dataset grid space. The mapping\n"
02383 " function determines how to map the surface values from many\n"
02384 " nodes to a single voxel.\n"
02385 "\n"
02386 " Surfaces (from the spec file) are specified using '-surf_A'\n"
02387 " (and '-surf_B', if a second surface is input). If two\n"
02388 " surfaces are input, then the computed segments over node\n"
02389 " pairs will be in the direction from surface A to surface B.\n"
02390 "\n"
02391 " The basic form of the algorithm is:\n"
02392 "\n"
02393 " o for each node pair (or single node)\n"
02394 " o form a segment based on the xyz node coordinates,\n"
02395 " adjusted by any '-f_pX_XX' options\n"
02396 " o divide the segment up into N steps, according to \n"
02397 " the '-f_steps' option\n"
02398 " o for each segment point\n"
02399 " o if the point is outside the space of the output\n"
02400 " dataset, skip it\n"
02401 " o locate the voxel in the output dataset which\n"
02402 " corresponds to this segment point\n"
02403 " o if the '-cmask' option was given, and the voxel\n"
02404 " is outside the implied mask, skip it\n"
02405 " o if the '-f_index' option is by voxel, and this\n"
02406 " voxel has already been considered, skip it\n"
02407 " o insert the surface node value, according to the\n"
02408 " user-specified '-map_func' option\n"
02409 "\n"
02410 " Surface Coordinates:\n"
02411 "\n"
02412 " Surface coordinates are assumed to be in the Dicom\n"
02413 " orientation. This information may come from the option\n"
02414 " pair of '-spec' and '-sv', with which the user provides\n"
02415 " the name of the SPEC FILE and the SURFACE VOLUME, along\n"
02416 " with '-surf_A' and optionally '-surf_B', used to specify\n"
02417 " actual surfaces by name. Alternatively, the surface\n"
02418 " coordinates may come from the '-surf_xyz_1D' option.\n"
02419 " See these option descriptions below.\n"
02420 "\n"
02421 " Note that the user must provide either the three options\n"
02422 " '-spec', '-sv' and '-surf_A', or the single option,\n"
02423 " '-surf_xyz_1D'.\n"
02424 "\n"
02425 " Surface Data:\n"
02426 "\n"
02427 " Surface domain data can be input via the '-sdata_1D'\n"
02428 " option. In such a case, the data is with respect to the\n"
02429 " input surface. The first column of the sdata_1D file\n"
02430 " should be a node index, and following columns are that\n"
02431 " node's data. See the '-sdata_1D' option for more info.\n"
02432 "\n"
02433 " If the surfaces have V values per node (pair), then the\n"
02434 " resulting AFNI dataset will have V sub-bricks (unless the\n"
02435 " user applies the '-data_expr' option).\n"
02436 "\n"
02437 " Mapping Functions:\n"
02438 "\n"
02439 " Mapping functions exist because a single volume voxel may\n"
02440 " be occupied by multiple surface nodes or segment points.\n"
02441 " Depending on how dense the surface mesh is, the number of\n"
02442 " steps provided by the '-f_steps' option, and the indexing\n"
02443 " type from '-f_index', even a voxel which is only 1 cubic\n"
02444 " mm in volume may have quite a few contributing points.\n"
02445 "\n"
02446 " The mapping function defines how multiple surface values\n"
02447 " are combined to get a single result in each voxel. For\n"
02448 " example, the 'max' function will take the maximum of all\n"
02449 " surface values contributing to each given voxel.\n"
02450 "\n"
02451 " Current mapping functions are listed under the '-map_func'\n"
02452 " option, below.\n"
02453 "\n"
02454 "------------------------------------------------------------\n"
02455 "\n"
02456 " examples:\n"
02457 "\n"
02458 " 1. Map a single surface to an anatomical volume domain,\n"
02459 " creating a simple mask of the surface. The output\n"
02460 " dataset will be fred_surf+orig, and the orientation and\n"
02461 " grid spacing will follow that of the grid parent. The\n"
02462 " output voxels will be 1 where the surface exists, and 0\n"
02463 " elsewhere.\n"
02464 "\n"
02465 " %s \\\n"
02466 " -spec fred.spec \\\n"
02467 " -surf_A pial \\\n"
02468 " -sv fred_anat+orig \\\n"
02469 " -grid_parent fred_anat+orig \\\n"
02470 " -map_func mask \\\n"
02471 " -prefix fred_surf\n"
02472 "\n"
02473 " 2. Map the cortical grey ribbon (between the white matter\n"
02474 " surface and the pial surface) to an AFNI volume, where\n"
02475 " the resulting volume is restriced to the mask implied by\n"
02476 " the -cmask option.\n"
02477 "\n"
02478 " Surface data will come from the file sdata_10.1D, which\n"
02479 " has 10 values per node, and lists only a portion of the\n"
02480 " entire set of surface nodes. Each node pair will be form\n"
02481 " a segment of 15 equally spaced points, the values from\n"
02482 " which will be applied to the output dataset according to\n"
02483 " the 'ave' filter. Since the index is over points, each\n"
02484 " of the 15 points will have its value applied to the\n"
02485 " appropriate voxel, even multiple times. This weights the\n"
02486 " resulting average by the fraction of each segment that\n"
02487 " occupies a given voxel.\n"
02488 "\n"
02489 " The output dataset will have 10 sub-bricks, according to\n"
02490 " the 10 values per node index in sdata_10.1D.\n"
02491 "\n"
02492 " %s \\\n"
02493 " -spec fred.spec \\\n"
02494 " -surf_A smoothwm \\\n"
02495 " -surf_B pial \\\n"
02496 " -sv fred_anat+orig \\\n"
02497 " -grid_parent 'fred_func+orig[0]' \\\n"
02498 " -cmask '-a fred_func+orig[2] -expr step(a-0.6)' \\\n"
02499 " -sdata_1D sdata_10.1D \\\n"
02500 " -map_func ave \\\n"
02501 " -f_steps 15 \\\n"
02502 " -f_index points \\\n"
02503 " -prefix fred_surf_ave\n"
02504 "\n"
02505 " 3. The inputs in this example are identical to those in\n"
02506 " example 2, including the surface dataset, sdata_10.1D.\n"
02507 " Again, the output dataset will have 10 sub-bricks.\n"
02508 "\n"
02509 " The surface values will be applied via the 'max_abs'\n"
02510 " filter, with the intention of assigning to each voxel the\n"
02511 " node value with the most significance. Here, the index\n"
02512 " method does not matter, so it is left as the default,\n"
02513 " 'voxel'.\n"
02514 "\n"
02515 " In this example, each node pair segment will be extended\n"
02516 " by 20%% into the white matter, and by 10%% outside of the\n"
02517 " grey matter, generating a \"thicker\" result.\n"
02518 "\n"
02519 " %s \\\n"
02520 " -spec fred.spec \\\n"
02521 " -surf_A smoothwm \\\n"
02522 " -surf_B pial \\\n"
02523 " -sv fred_anat+orig \\\n"
02524 " -grid_parent 'fred_func+orig[0]' \\\n"
02525 " -cmask '-a fred_func+orig[2] -expr step(a-0.6)' \\\n"
02526 " -sdata_1D sdata_10.1D \\\n"
02527 " -map_func max_abs \\\n"
02528 " -f_steps 15 \\\n"
02529 " -f_p1_fr -0.2 \\\n"
02530 " -f_pn_fr 0.1 \\\n"
02531 " -prefix fred_surf_max_abs\n"
02532 "\n"
02533 " 4. This is simliar to example 2. Here, the surface nodes\n"
02534 " (coordinates) come from 'surf_coords_2.1D'. But these\n"
02535 " coordinates do not happen to be in Dicomm orientation,\n"
02536 " they are in the same orientation as the grid parent, so\n"
02537 " the '-sxyz_orient_as_gpar' option is applied.\n"
02538 "\n"
02539 " Even though the data comes from 'sdata_10.1D', the output\n"
02540 " AFNI dataset will only have 1 sub-brick. That is because\n"
02541 " of the '-data_expr' option. Here, each applied surface\n"
02542 " value will be the average of the sines of the first 3\n"
02543 " data values (columns of sdata_10.1D).\n"
02544 "\n"
02545 " %s \\\n"
02546 " -surf_xyz_1D surf_coords_2.1D \\\n"
02547 " -sxyz_orient_as_gpar \\\n"
02548 " -grid_parent 'fred_func+orig[0]' \\\n"
02549 " -sdata_1D sdata_10.1D \\\n"
02550 " -data_expr '(sin(a)+sin(b)+sin(c))/3' \\\n"
02551 " -map_func ave \\\n"
02552 " -f_steps 15 \\\n"
02553 " -f_index points \\\n"
02554 " -prefix fred_surf_ave_sine\n"
02555 "\n"
02556 " 5. In this example, voxels will get the maximum value from\n"
02557 " column 3 of sdata_10.1D (as usual, column 0 is used for\n"
02558 " node indices). The output dataset will have 1 sub-brick.\n"
02559 "\n"
02560 " Here, the output dataset is forced to be of type 'short',\n"
02561 " regardless of what the grid parent is. Also, there will\n"
02562 " be no scaling factor applied.\n"
02563 "\n"
02564 " To track the numbers for surface node #1234, the '-dnode'\n"
02565 " option has been used, along with '-debug'. Additionally,\n"
02566 " '-dvoxel' is used to track the results for voxel #6789.\n"
02567 "\n"
02568 " %s \\\n"
02569 " -spec fred.spec \\\n"
02570 " -surf_A smoothwm \\\n"
02571 " -surf_B pial \\\n"
02572 " -sv fred_anat+orig \\\n"
02573 " -grid_parent 'fred_func+orig[0]' \\\n"
02574 " -sdata_1D sdata_10.1D'[0,3]' \\\n"
02575 " -map_func max \\\n"
02576 " -f_steps 15 \\\n"
02577 " -datum short \\\n"
02578 " -noscale \\\n"
02579 " -debug 2 \\\n"
02580 " -dnode 1234 \\\n"
02581 " -dvoxel 6789 \\\n"
02582 " -prefix fred_surf_max\n"
02583 "\n"
02584 "------------------------------------------------------------\n"
02585 "\n"
02586 " REQUIRED COMMAND ARGUMENTS:\n"
02587 "\n"
02588 " -spec SPEC_FILE : SUMA spec file\n"
02589 "\n"
02590 " e.g. -spec fred.spec\n"
02591 "\n"
02592 " The surface specification file contains the list of\n"
02593 " mappable surfaces that are used.\n"
02594 "\n"
02595 " See @SUMA_Make_Spec_FS and @SUMA_Make_Spec_SF.\n"
02596 "\n"
02597 " Note: this option, along with '-sv', may be replaced\n"
02598 " by the '-surf_xyz_1D' option.\n"
02599 "\n"
02600 " -surf_A SURF_NAME : specify surface A (from spec file)\n"
02601 " -surf_B SURF_NAME : specify surface B (from spec file)\n"
02602 "\n"
02603 " e.g. -surf_A smoothwm\n"
02604 " e.g. -surf_A lh.smoothwm\n"
02605 " e.g. -surf_B lh.pial\n"
02606 "\n"
02607 " This parameter is used to tell the program with surfaces\n"
02608 " to use. The '-surf_A' parameter is required, but the\n"
02609 " '-surf_B' parameter is an option.\n"
02610 "\n"
02611 " The surface names must uniquely match those in the spec\n"
02612 " file, though a sub-string match is good enough. The\n"
02613 " surface names are compared with the names of the surface\n"
02614 " node coordinate files.\n"
02615 "\n"
02616 " For instance, given a spec file that has only the left\n"
02617 " hemisphere in it, 'pial' should produce a unique match\n"
02618 " with lh.pial.asc. But if both hemispheres are included,\n"
02619 " then 'pial' would not be unique (matching rh.pail.asc,\n"
02620 " also). In that case, 'lh.pial' would be better.\n"
02621 "\n"
02622 " -sv SURFACE_VOLUME : AFNI dataset\n"
02623 "\n"
02624 " e.g. -sv fred_anat+orig\n"
02625 "\n"
02626 " This is the AFNI dataset that the surface is mapped to.\n"
02627 " This dataset is used for the initial surface node to xyz\n"
02628 " coordinate mapping, in the Dicom orientation.\n"
02629 "\n"
02630 " Note: this option, along with '-spec', may be replaced\n"
02631 " by the '-surf_xyz_1D' option.\n"
02632 "\n"
02633 " -surf_xyz_1D SXYZ_NODE_FILE : 1D coordinate file\n"
02634 "\n"
02635 " e.g. -surf_xyz_1D my_surf_coords.1D\n"
02636 "\n"
02637 " This ascii file contains a list of xyz coordinates to be\n"
02638 " considered as a surface, or 2 sets of xyz coordinates to\n"
02639 " considered as a surface pair. As usual, these points\n"
02640 " are assumed to be in Dicom orientation. Another option\n"
02641 " for coordinate orientation is to use that of the grid\n"
02642 " parent dataset. See '-sxyz_orient_as_gpar' for details.\n"
02643 "\n"
02644 " This option is an alternative to the pair of options, \n"
02645 " '-spec' and '-sv'.\n"
02646 "\n"
02647 " The number of rows of the file should equal the number\n"
02648 " of nodes on each surface. The number of columns should\n"
02649 " be either 3 for a single surface, or 6 for two surfaces.\n"
02650 " \n"
02651 " sample line of an input file (one surface):\n"
02652 " \n"
02653 " 11.970287 2.850751 90.896111\n"
02654 " \n"
02655 " sample line of an input file (two surfaces):\n"
02656 " \n"
02657 " 11.97 2.85 90.90 12.97 2.63 91.45\n"
02658 " \n"
02659 "\n"
02660 " -grid_parent AFNI_DSET : AFNI dataset\n"
02661 "\n"
02662 " e.g. -grid_parent fred_function+orig\n"
02663 "\n"
02664 " This dataset is used as a grid and orientation master\n"
02665 " for the output AFNI dataset.\n"
02666 "\n"
02667 " -map_func MAP_FUNC : surface to dataset function\n"
02668 "\n"
02669 " e.g. -map_func max\n"
02670 " e.g. -map_func mask -f_steps 20\n"
02671 "\n"
02672 " This function applies to the case where multiple data\n"
02673 " points get mapped to a single voxel, which is expected\n"
02674 " since surfaces tend to have a much higher resolution\n"
02675 " than AFNI volumes. In the general case data points come\n"
02676 " from each point on each partitioned line segment, with\n"
02677 " one segment per node pair. Note that these segments may\n"
02678 " have length zero, such as when only a single surface is\n"
02679 " input.\n"
02680 "\n"
02681 " See \"Mapping Functions\" above, for more information.\n"
02682 "\n"
02683 " The current mapping function for one surface is:\n"
02684 "\n"
02685 " mask : For each xyz location, set the corresponding\n"
02686 " voxel to 1.\n"
02687 "\n"
02688 " The current mapping functions for two surfaces are as\n"
02689 " follows. These descriptions are per output voxel, and\n"
02690 " over the values of all points mapped to a given voxel.\n"
02691 "\n"
02692 " mask2 : if any points are mapped to the voxel, set\n"
02693 " the voxel value to 1\n"
02694 "\n"
02695 " ave : average all values\n"
02696 "\n"
02697 " count : count the number of mapped data points\n"
02698 "\n"
02699 " min : find the minimum value from all mapped points\n"
02700 "\n"
02701 " max : find the maximum value from all mapped points\n"
02702 "\n"
02703 " max_abs: find the number with maximum absolute value\n"
02704 " (the resulting value will retain its sign)\n"
02705 "\n"
02706 " -prefix OUTPUT_PREFIX : prefix for the output dataset\n"
02707 "\n"
02708 " e.g. -prefix anat_surf_mask\n"
02709 "\n"
02710 " This is used to specify the prefix of the resulting AFNI\n"
02711 " dataset.\n"
02712 "\n"
02713 " ------------------------------\n"
02714 " SUB-SURFACE DATA FILE OPTIONS:\n"
02715 "\n"
02716 " -sdata_1D SURF_DATA.1D : 1D sub-surface file, with data\n"
02717 "\n"
02718 " e.g. -sdata_1D roi3.1D\n"
02719 "\n"
02720 " This is used to specify a 1D file, which contains\n"
02721 " surface indices and data. The indices refer to the\n"
02722 " surface(s) read from the spec file.\n"
02723 " \n"
02724 " The format of this data file is a surface index and a\n"
02725 " list of data values on each row. To be a valid 1D file,\n"
02726 " each row must have the same number of columns.\n"
02727 "\n"
02728 " ------------------------------\n"
02729 " OPTIONS SPECIFIC TO SEGMENT SELECTION:\n"
02730 "\n"
02731 " (see \"The basic form of the algorithm\" for more details)\n"
02732 "\n"
02733 " -f_steps NUM_STEPS : partition segments\n"
02734 "\n"
02735 " e.g. -f_steps 10\n"
02736 " default: -f_steps 2 (or 1, the number of surfaces)\n"
02737 "\n"
02738 " This option specifies the number of points to divide\n"
02739 " each line segment into, before mapping the points to the\n"
02740 " AFNI volume domain. The default is the number of input\n"
02741 " surfaces (usually, 2). The default operation is to have\n"
02742 " the segment endpoints be the actual surface nodes,\n"
02743 " unless they are altered with the -f_pX_XX options.\n"
02744 "\n"
02745 " -f_index TYPE : index by points or voxels\n"
02746 "\n"
02747 " e.g. -f_index points\n"
02748 " e.g. -f_index voxels\n"
02749 " default: -f_index voxels\n"
02750 "\n"
02751 " Along a single segment, the default operation is to\n"
02752 " apply only those points mapping to a new voxel. The\n"
02753 " effect of the default is that a given voxel will have\n"
02754 " at most one value applied per voxel pair.\n"
02755 "\n"
02756 " If the user applies this option with 'points' or 'nodes'\n"
02757 " as the argument, then every point along the segment will\n"
02758 " be applied. This may be preferred if, for example, the\n"
02759 " user wishes to have the average weighted by the number\n"
02760 " of points occupying a voxel, not just the number of node\n"
02761 " pair segments.\n"
02762 "\n"
02763 " Note: the following -f_pX_XX options are used to alter the\n"
02764 " locations of the segment endpoints, per node pair.\n"
02765 " The segments are directed, from the node on the first\n"
02766 " surface to the node on the second surface. To modify\n"
02767 " the first endpoint, use a -f_p1_XX option, and use\n"
02768 " -f_pn_XX to modify the second.\n"
02769 "\n"
02770 " -f_p1_fr FRACTION : offset p1 by a length fraction\n"
02771 "\n"
02772 " e.g. -f_p1_fr -0.2\n"
02773 " e.g. -f_p1_fr -0.2 -f_pn_fr 0.2\n"
02774 "\n"
02775 " This option moves the first endpoint, p1, by a distance\n"
02776 " of the FRACTION times the original segment length. If\n"
02777 " the FRACTION is positive, it moves in the direction of\n"
02778 " the second endpoint, pn.\n"
02779 "\n"
02780 " In the example, p1 is moved by 20%% away from pn, which\n"
02781 " will increase the length of each segment.\n"
02782 "\n"
02783 " -f_pn_fr FRACTION : offset pn by a length fraction\n"
02784 "\n"
02785 " e.g. -f_pn_fr 0.2\n"
02786 " e.g. -f_p1_fr -0.2 -f_pn_fr 0.2\n"
02787 "\n"
02788 " This option moves pn by a distance of the FRACTION times\n"
02789 " the original segment length, in the direction from p1 to\n"
02790 " pn. So a positive fraction extends the segment, and a\n"
02791 " negative fraction reduces it.\n"
02792 "\n"
02793 " In the example above, using 0.2 adds 20%% to the segment\n"
02794 " length past the original pn.\n"
02795 "\n"
02796 " -f_p1_mm DISTANCE : offset p1 by a distance in mm.\n"
02797 "\n"
02798 " e.g. -f_p1_mm -1.0\n"
02799 " e.g. -f_p1_mm -1.0 -f_pn_fr 1.0\n"
02800 "\n"
02801 " This option moves p1 by DISTANCE mm., in the direction\n"
02802 " of pn. If the DISTANCE is positive, the segment gets\n"
02803 " shorter. If DISTANCE is negative, the segment will get\n"
02804 " longer.\n"
02805 "\n"
02806 " In the example, p1 is moved away from pn, extending the\n"
02807 " segment by 1 millimeter.\n"
02808 "\n"
02809 " -f_pn_mm DISTANCE : offset pn by a distance in mm.\n"
02810 "\n"
02811 " e.g. -f_pn_mm 1.0\n"
02812 " e.g. -f_p1_mm -1.0 -f_pn_fr 1.0\n"
02813 "\n"
02814 " This option moves pn by DISTANCE mm., in the direction\n"
02815 " from the first point to the second. So if DISTANCE is\n"
02816 " positive, the segment will get longer. If DISTANCE is\n"
02817 " negative, the segment will get shorter.\n"
02818 "\n"
02819 " In the example, pn is moved 1 millimeter farther from\n"
02820 " p1, extending the segment by that distance.\n"
02821 "\n"
02822 " ------------------------------\n"
02823 " GENERAL OPTIONS:\n"
02824 "\n"
02825 " -cmask MASK_COMMAND : command for dataset mask\n"
02826 "\n"
02827 " e.g. -cmask '-a fred_func+orig[2] -expr step(a-0.8)'\n"
02828 "\n"
02829 " This option will produce a mask to be applied to the\n"
02830 " output dataset. Note that this mask should form a\n"
02831 " single sub-brick.\n"
02832 "\n"
02833 " This option follows the style of 3dmaskdump (since the\n"
02834 " code for it was, uh, borrowed from there (thanks Bob!)).\n"
02835 "\n"
02836 " See '3dmaskdump -help' for more information.\n"
02837 "\n"
02838 " -data_expr EXPRESSION : apply expression to surface input\n"
02839 "\n"
02840 " e.g. -data_expr 17\n"
02841 " e.g. -data_expr '(a+b+c+d)/4'\n"
02842 " e.g. -data_expr '(sin(a)+sin(b))/2'\n"
02843 "\n"
02844 " This expression is applied to the list of data values\n"
02845 " from the surface data file input via '-sdata_1D'. The\n"
02846 " expression is applied for each node or node pair, to the\n"
02847 " list of data values corresponding to that node.\n"
02848 "\n"
02849 " The letters 'a' through 'z' may be used as input, and\n"
02850 " refer to columns 1 through 26 of the data file (where\n"
02851 " column 0 is a surface node index). The data file must\n"
02852 " have enough columns to support the expression. Is is\n"
02853 " valid to have a constant expression without a data file.\n"
02854 "\n"
02855 " -datum DTYPE : set data type in output dataset\n"
02856 "\n"
02857 " e.g. -datum short\n"
02858 " default: same as that of grid parent\n"
02859 "\n"
02860 " This option specifies the data type for the output AFNI\n"
02861 " dataset. Valid choices are byte, short and float, which\n"
02862 " are 1, 2 and 4 bytes for each data point, respectively.\n"
02863 "\n"
02864 " -debug LEVEL : verbose output\n"
02865 "\n"
02866 " e.g. -debug 2\n"
02867 "\n"
02868 " This option is used to print out status information \n"
02869 " during the execution of the program. Current levels are\n"
02870 " from 0 to 5.\n"
02871 "\n"
02872 " -dnode DEBUG_NODE : extra output for that node\n"
02873 "\n"
02874 " e.g. -dnode 123456\n"
02875 "\n"
02876 " This option requests additional debug output for the\n"
02877 " given surface node. This index is with respect to the\n"
02878 " input surface (included in the spec file, or through the\n"
02879 " '-surf_xyz_1D' option).\n"
02880 "\n"
02881 " This will have no effect without the '-debug' option.\n"
02882 "\n"
02883 " -dvoxel DEBUG_VOXEL : extra output for that voxel\n"
02884 "\n"
02885 " e.g. -dvoxel 234567\n"
02886 "\n"
02887 " This option requests additional debug output for the\n"
02888 " given volume voxel. This 1-D index is with respect to\n"
02889 " the output AFNI dataset. One good way to find a voxel\n"
02890 " index to supply is from output via the '-dnode' option.\n"
02891 "\n"
02892 " This will have no effect without the '-debug' option.\n"
02893 "\n"
02894 " -hist : show revision history\n"
02895 "\n"
02896 " Display module history over time.\n"
02897 "\n"
02898 " -help : show this help\n"
02899 "\n"
02900 " If you can't get help here, please get help somewhere.\n"
02901 "\n"
02902 " -noscale : no scale factor in output dataset\n"
02903 "\n"
02904 " If the output dataset is an integer type (byte, shorts\n"
02905 " or ints), then the output dataset may end up with a\n"
02906 " scale factor attached (see 3dcalc -help). With this\n"
02907 " option, the output dataset will not be scaled.\n"
02908 "\n"
02909 " -sxyz_orient_as_gpar : assume gpar orientation for sxyz\n"
02910 "\n"
02911 " This option specifies that the surface coordinate points\n"
02912 " in the '-surf_xyz_1D' option file have the orientation\n"
02913 " of the grid parent dataset.\n"
02914 "\n"
02915 " When the '-surf_xyz_1D' option is applied the surface\n"
02916 " coordinates are assumed to be in Dicom orientation, by\n"
02917 " default. This '-sxyz_orient_as_gpar' option overrides\n"
02918 " the Dicom default, specifying that the node coordinates\n"
02919 " are in the same orientation as the grid parent dataset.\n"
02920 "\n"
02921 " See the '-surf_xyz_1D' option for more information.\n"
02922 "\n"
02923 " -version : show version information\n"
02924 "\n"
02925 " Show version and compile date.\n"
02926 "\n"
02927 "------------------------------------------------------------\n"
02928 "\n"
02929 " Author: R. Reynolds - %s\n"
02930 "\n"
02931 " (many thanks to Z. Saad and R.W. Cox)\n"
02932 "\n",
02933 prog, prog,
02934 prog, prog, prog, prog, prog,
02935 VERSION );
02936 }
02937 else if ( level == S2V_USE_HIST )
02938 fputs(g_history, stdout);
02939 else if ( level == S2V_USE_VERSION )
02940 printf( "%s : %s, compile date: %s\n", prog, VERSION, __DATE__ );
02941 else
02942 fprintf( stderr, "usage called with illegal level <%d>\n", level );
02943
02944 RETURN(-1);
02945 }
02946
02947
02948
02949
02950
02951
02952 int disp_opts_t ( char * info, opts_t * opts )
02953 {
02954 ENTRY("disp_opts_t");
02955
02956 if ( info )
02957 fputs( info, stderr );
02958
02959 if ( opts == NULL )
02960 {
02961 fprintf(stderr, "disp_opts_t: opts == NULL\n" );
02962 RETURN(-1);
02963 }
02964
02965 fprintf(stderr,
02966 "options struct at %p :\n"
02967 " gpar_file = %s\n"
02968 " oset_file = %s\n"
02969 " spec_file = %s\n"
02970 " sv_file = %s\n"
02971 " surf_xyz_1D_file = %s\n"
02972 " sdata_file_1D = %s\n"
02973 " sdata_file_niml = %s\n"
02974 " cmask_cmd = %s\n"
02975 " data_expr = %s\n"
02976 " map_str = %s\n"
02977 " datum_str = %s\n"
02978 " f_index_str = %s\n"
02979 " sxyz_ori_gpar = %d\n"
02980 " debug, dnode, dvox = %d, %d, %d\n"
02981 " noscale, f_steps = %d, %d\n"
02982 " f_p1_fr, f_pn_fr = %f, %f\n"
02983 " f_p1_mm, f_pn_mm = %f, %f\n"
02984 , opts,
02985 CHECK_NULL_STR(opts->gpar_file), CHECK_NULL_STR(opts->oset_file),
02986 CHECK_NULL_STR(opts->spec_file), CHECK_NULL_STR(opts->sv_file),
02987 CHECK_NULL_STR(opts->surf_xyz_1D_file),
02988 CHECK_NULL_STR(opts->sdata_file_1D),
02989 CHECK_NULL_STR(opts->sdata_file_niml),
02990 CHECK_NULL_STR(opts->cmask_cmd), CHECK_NULL_STR(opts->data_expr),
02991 CHECK_NULL_STR(opts->map_str), CHECK_NULL_STR(opts->datum_str),
02992 CHECK_NULL_STR(opts->f_index_str), opts->sxyz_ori_gpar,
02993 opts->debug, opts->dnode, opts->dvox, opts->noscale, opts->f_steps,
02994 opts->f_p1_fr, opts->f_pn_fr, opts->f_p1_mm, opts->f_pn_mm
02995 );
02996
02997 RETURN(0);
02998 }
02999
03000
03001
03002
03003
03004
03005 int disp_param_t ( char * info, param_t * p )
03006 {
03007 ENTRY("disp_param_t");
03008
03009 if ( info )
03010 fputs( info, stderr );
03011
03012 if ( p == NULL )
03013 {
03014 fprintf(stderr,"disp_param_t: p == NULL\n");
03015 RETURN(-1);
03016 }
03017
03018 fprintf(stderr,
03019 "param_t struct at %p :\n"
03020 " gpar : vcheck = %p : %s\n"
03021 " oset : vcheck = %p : %s\n"
03022 " sxyz_im, sdata_im = %p, %p\n"
03023 " f3mm_min (xyz) = (%f, %f, %f)\n"
03024 " f3mm_max (xyz) = (%f, %f, %f)\n"
03025 " nvox, nsubs = %d, %d\n"
03026 " cmask = %p\n"
03027 " ncmask, ccount = %d, %d\n"
03028 , p,
03029 p->gpar, ISVALID_DSET(p->gpar) ? "valid" : "invalid",
03030 p->oset, ISVALID_DSET(p->oset) ? "valid" : "invalid",
03031 p->sxyz_im, p->sdata_im,
03032 p->f3mm_min.xyz[0], p->f3mm_min.xyz[1], p->f3mm_min.xyz[2],
03033 p->f3mm_max.xyz[0], p->f3mm_max.xyz[1], p->f3mm_max.xyz[2],
03034 p->nvox, p->nsubs, p->cmask, p->ncmask, p->ccount
03035 );
03036
03037 RETURN(0);
03038 }
03039
03040
03041
03042
03043
03044
03045 int disp_node_list_t ( char * info, node_list_t * d )
03046 {
03047 ENTRY("disp_node_list_t");
03048
03049 if ( info )
03050 fputs( info, stderr );
03051
03052 if ( d == NULL )
03053 {
03054 fprintf(stderr,"disp_node_list_t: d == NULL\n");
03055 RETURN(-1);
03056 }
03057
03058 fprintf(stderr,
03059 "node_list_t struct at %p :\n"
03060 " nodes = %p\n"
03061 " depth, nnodes = %d, %d\n"
03062 " fdata, ilist = %p, %p\n"
03063 " ilen = %d\n"
03064 , d,
03065 d->nodes, d->depth, d->nnodes,
03066 d->fdata, d->ilist, d->ilen
03067 );
03068
03069 RETURN(0);
03070 }
03071
03072
03073
03074
03075
03076
03077 int disp_s2v_opts_t ( char * info, s2v_opts_t * sopt )
03078 {
03079 ENTRY("disp_s2v_opts_t");
03080
03081 if ( info )
03082 fputs( info, stderr );
03083
03084 if ( sopt == NULL )
03085 {
03086 fprintf(stderr,"disp_s2v_opts_t: sopt == NULL\n" );
03087 RETURN(-1);
03088 }
03089
03090 fprintf(stderr,
03091 "s2v_opts_t struct at %p :\n"
03092 " map, datum, noscale = %d, %d, %d\n"
03093 " debug, dnode, dvox = %d, %d, %d\n"
03094 " sxyz_ori_gpar = %d\n"
03095 " cmask = %p\n"
03096 " f_steps, f_index = %d, %d\n"
03097 " f_p1_fr, f_pn_fr = %f, %f\n"
03098 " f_p1_mm, f_pn_mm = %f, %f\n"
03099 , sopt,
03100 sopt->map, sopt->datum, sopt->noscale,
03101 sopt->debug, sopt->dnode, sopt->dvox, sopt->sxyz_ori_gpar,
03102 sopt->cmask, sopt->f_steps, sopt->f_index,
03103 sopt->f_p1_fr, sopt->f_pn_fr, sopt->f_p1_mm, sopt->f_pn_mm
03104 );
03105
03106 RETURN(0);
03107 }
03108
03109
03110
03111
03112
03113
03114 int disp_parser_t ( char * info, parser_t * d )
03115 {
03116 int c;
03117
03118 ENTRY("disp_parser_t");
03119
03120 if ( info )
03121 fputs( info, stderr );
03122
03123 if ( d == NULL )
03124 {
03125 fprintf(stderr,"disp_parser_t: d == NULL\n" );
03126 RETURN(-1);
03127 }
03128
03129 fprintf(stderr, "parser_t struct at %p :\n"
03130 " pcode = %p\n"
03131 " max_sym = %d\n",
03132 d, d->pcode, d->max_sym );
03133
03134 if ( d->pcode )
03135 {
03136 fprintf(stderr, " num_code = %d\n"
03137 " c_code = %s\n",
03138 d->pcode->num_code, d->pcode->c_code );
03139
03140 fprintf(stderr, " has_sym =");
03141 for ( c = 0; c < 26; c++ )
03142 fprintf(stderr, " %d", d->has_sym[c]);
03143 fputc('\n', stderr);
03144 }
03145
03146 RETURN(0);
03147 }
03148
03149
03150
03151
03152
03153
03154
03155
03156 int set_3dmm_bounds ( THD_3dim_dataset *dset, THD_fvec3 *min, THD_fvec3 *max)
03157 {
03158 float tmp;
03159 int c;
03160
03161 ENTRY("set_3dmm_bounds");
03162
03163 if ( !dset || !min || !max )
03164 {
03165 fprintf(stderr, "** invalid params to set_3dmm_bounds: (%p,%p,%p)\n",
03166 dset, min, max );
03167 RETURN(-1);
03168 }
03169
03170
03171 min->xyz[0] = DSET_XORG(dset) - 0.5 * DSET_DX(dset);
03172 max->xyz[0] = min->xyz[0] + DSET_NX(dset) * DSET_DX(dset);
03173
03174 min->xyz[1] = DSET_YORG(dset) - 0.5 * DSET_DY(dset);
03175 max->xyz[1] = min->xyz[1] + DSET_NY(dset) * DSET_DY(dset);
03176
03177 min->xyz[2] = DSET_ZORG(dset) - 0.5 * DSET_DZ(dset);
03178 max->xyz[2] = min->xyz[2] + DSET_NZ(dset) * DSET_DZ(dset);
03179
03180 for ( c = 0; c < 3; c++ )
03181 if ( min->xyz[c] > max->xyz[c] )
03182 {
03183 tmp = min->xyz[c];
03184 min->xyz[c] = max->xyz[c];
03185 max->xyz[c] = tmp;
03186 }
03187
03188 RETURN(0);
03189 }
03190
03191
03192
03193
03194
03195
03196 float dist_f3mm( THD_fvec3 * p1, THD_fvec3 * p2 )
03197 {
03198 double d0, d1, d2;
03199
03200 ENTRY("dist_f3mm");
03201
03202 if ( p1 == NULL || p2 == NULL )
03203 {
03204 fprintf( stderr, "** dist_f3mm: invalid params (%p,%p)\n", p1, p2 );
03205 RETURN(0.0);
03206 }
03207
03208 d0 = p1->xyz[0] - p2->xyz[0];
03209 d1 = p1->xyz[1] - p2->xyz[1];
03210 d2 = p1->xyz[2] - p2->xyz[2];
03211
03212 RETURN(sqrt(d0*d0 + d1*d1 + d2*d2));
03213 }
03214
03215
03216
03217
03218
03219
03220
03221 int f3mm_out_of_bounds( THD_fvec3 * cp, THD_fvec3 * min, THD_fvec3 * max )
03222 {
03223 int c;
03224
03225 ENTRY("f3mm_out_of_bounds");
03226
03227 if ( !cp || !min || !max )
03228 RETURN(-1);
03229
03230 for ( c = 0; c < 3; c++ )
03231 {
03232 if ( ( cp->xyz[c] < min->xyz[c] ) ||
03233 ( cp->xyz[c] > max->xyz[c] ) )
03234 RETURN(-1);
03235 }
03236
03237 RETURN(0);
03238 }
03239
03240
03241
03242
03243
03244
03245
03246 int integral_doubles( double * dp, int nvals )
03247 {
03248 ENTRY("integral_doubles");
03249
03250 while ( nvals > 0 )
03251 {
03252 if ( ((double)(int)*dp) != *dp )
03253 RETURN(0);
03254
03255 dp++;
03256 nvals--;
03257 }
03258
03259 RETURN(1);
03260 }
03261
