/*----------------------------------------------------------------------
* 3dSurf2Vol - create an AFNI volume dataset from a surface
*
* This program is meant to take as input a surface and an AFNI dataset,
* and to output a new AFNI dataset consisting of the surface mapped to
* the dataset grid space.
*
* The surface points are mapped to xyz coordinates, according to the
* SURF_VOL (surface volume) AFNI dataset. These coordinates are then
* matched to voxels in the grid parent AFNI dataset. The output dataset
* will have the same grid as the grid parent dataset, with values coming
* from the surface, and subject to any mask (see -cmask).
*
* usage:
* 3dSurf2Vol [options] -spec SPEC_FILE -sv SURF_VOL -grid_parent AFNI_DSET
* -map_func MAP_FUNC -prefix OUTPUT_DSET
*
* options:
*
* -help
* -version
*
* -spec SPEC_FILE
* -surf_A SURF_NAME_A
* -surf_B SURF_NAME_B
* -sv SURF_VOL
* -grid_parent AFNI_DSET
* -prefix OUTPUT_DSET
* -map_func MAP_FUNC
*
* -surf_xyz_1D SURF_XYZ.1D
* -sdata_1D SURF_DATA.1D
*
* -cmask MASK_COMMAND
* -data_expr EXPRESSION
* -datum D_TYPE
* -debug LEVEL
* -dnode NODE_NUM
* -dvoxel VOXEL_NUM
* -noscale
*
* -f_steps SEGMENT_STEPS
* -f_index STEP_TYPE
* -f_p1_fr FRACTION
* -f_pn_fr FRACTION
* -f_p1_mm DISTANCE
* -f_pn_mm DISTANCE
*
* examples:
*
* 3dSurf2Vol -spec SubjectA.spec \
* -sv SubjectA_spgr+orig \
* -grid_parent SubjA_EPI+orig \
* -map_func mask2 \
* -prefix SubjA_mask
*
* 3dSurf2Vol -spec SubjectA.spec \
* -sv SubjectA_spgr+orig \
* -cmask '-a SubjA.func+orig[2] -expr step(a-0.6)' \
* -grid_parent SubjA_EPI+orig \
* -map_func mask2 \
* -debug 2 \
* -prefix SubjA_mask
*
*----------------------------------------------------------------------
*/
static char g_history[] =
"----------------------------------------------------------------------\n"
"history:\n"
"\n"
"1.0 May 29, 2003 [rickr]\n"
" - initial release\n"
"\n"
"1.1 June 11, 2003 [rickr]\n"
" - small reorg of s2v_fill_mask2() (should have no effect)\n"
" - improve description of -f_steps option\n"
"\n"
"1.2 July 21, 2003 [rickr]\n"
" - make sure input points fit in output dataset\n"
" - add min/max distance output, along with out-of-bounds count\n"
"\n"
"2.0 October 2, 2003 [rickr]\n"
" - Major changes accepting surface data, surface coordinates, output data\n"
" type, debug options, multiple sub-brick output, and node pair segment\n"
" alterations.\n"
" - Added many options: -surf_xyz_1D, -sdata_1D, -data_expr, -datum,\n"
" -dnode, -dvoxel, -f_index, -f_p1_fr, -f_pn_fr,\n"
" -f_p1_mm, -f_pn_mm\n"
"\n"
"2.1 October 20, 2003 [rickr]\n"
" - call the new engine function, SUMA_LoadSpec_eng()\n"
" (this will restrict the debug output from SUMA_LoadSpec())\n"
"\n"
"2.2 December 15, 2003 [rickr]\n"
" - added program arguments '-surf_A' and '-surf_B' (-surf_A is required)\n"
" - added option '-hist' (for program history)\n"
" - explicit pointer init to NULL (a.o.t. memset() to 0)\n"
"\n"
"3.0 December 18, 2003 [rickr]\n"
" - removed requirement of 2 surfaces for most functions\n"
" (so now all functions work with either 1 or 2 input surfaces)\n"
"\n"
"3.1 January 23, 2004 [rickr]\n"
" - SUMA_isINHmappable() is depricated, check with AnatCorrect field\n"
" - reversed order of output for '-hist' option\n"
"\n"
"3.2 February 10, 2004 [rickr]\n"
" - output a little more debug info for !AnatCorrect case\n"
"\n"
"3.3 March 26, 2004 [ziad]\n"
" - DsetList added to SUMA_LoadSpec_eng() call\n"
"\n"
"3.4 June 21, 2004 [rickr]\n"
" - Fixed -surf_xyz_1D option (broken in v3.0 on nsurf test).\n"
"\n"
"3.5 July 22, 2004 [rickr]\n"
" - Fixed bug with test for valid sdata_1D file.\n"
"\n"
"3.6 July 28, 2004 [rickr]\n"
" - Fixed bug: old change caused the default f_steps to revert to 1,\n"
" now it is set back to 2 (thanks, Kuba).\n"
"\n"
"3.6a March 22, 2005 [rickr]\n"
" - removed tabs\n"
"----------------------------------------------------------------------\n";
#define VERSION "version 3.6a (March 22, 2005)"
/*----------------------------------------------------------------------
* todo:
* - handle niml input
*----------------------------------------------------------------------
*/
#include "mrilib.h"
#include "parser.h"
#include "SUMA_suma.h"
#include "SUMA_3dSurf2Vol.h"
/* globals */
SUMA_SurfaceViewer * SUMAg_SVv = NULL; /* array of Surf View structs */
int SUMAg_N_SVv = 0; /* length of SVv array */
SUMA_DO * SUMAg_DOv = NULL; /* array of Displayable Objects */
int SUMAg_N_DOv = 0; /* length of DOv array */
SUMA_CommonFields * SUMAg_CF = NULL; /* info common to all viewers */
/* this must match s2v_map_num enum */
char * gs2v_map_names[] = { "none", "mask", "mask2", "ave", "count",
"min", "max", "max_abs" };
/* AFNI prototype */
extern void machdep( void );
#define MAIN
/*----------------------------------------------------------------------*/
int main( int argc , char * argv[] )
{
SUMA_SurfSpecFile spec;
node_list_t node_list;
param_t params;
s2v_opts_t sopt;
opts_t opts;
int ret_val;
mainENTRY("3dSurf2Vol main");
machdep();
AFNI_logger("3dSurf2Vol",argc,argv);
/* validate inputs and init options structure */
if ( ( ret_val = init_options(&opts, argc, argv) ) != 0 )
return ret_val;
if ( ( ret_val = validate_options(&opts, ¶ms) ) != 0 )
return ret_val;
if ( (ret_val = set_smap_opts( &opts, ¶ms, &sopt )) != 0 )
return ret_val;
/* initialize the spec ZSS Jan 9 06*/
if (!SUMA_AllocSpecFields(&spec)) {
fprintf( stderr, "** failed to initialize spec\n" );
return(-1);
}
/* read surface files */
ret_val = read_surf_files(&opts, ¶ms, &spec, &sopt, &node_list);
if ( ret_val == 0 )
ret_val = write_output( &sopt, &opts, ¶ms, &node_list, argc, argv );
/* free memory */
final_clean_up( &node_list );
/* free the spec ZSS Jan 9 06*/
if (!SUMA_FreeSpecFields(&spec)) {
fprintf( stderr, "** failed to free spec\n" );
return(-1);
}
return ret_val;
}
/*----------------------------------------------------------------------
* write_output - create and write a new afni dataset
*----------------------------------------------------------------------
*/
int write_output ( s2v_opts_t * sopt, opts_t * opts, param_t * p,
node_list_t * N, int argc, char * argv[] )
{
ENTRY("write_output");
if ( sopt == NULL || opts == NULL || p == NULL || N == NULL )
{
fprintf( stderr, "** s2v_wo - bad params (%p,%p,%p,%p)\n",
sopt, opts, p, N );
RETURN(-1);
}
p->oset = s2v_nodes2volume( N, p, sopt );
if ( p->oset == NULL )
RETURN(-1);
EDIT_dset_items( p->oset, ADN_prefix, opts->oset_file, ADN_none );
if ( THD_is_file(DSET_HEADNAME(p->oset)) )
{
fprintf( stderr, "** cannot overwrite existing dataset '%s'\n",
DSET_HEADNAME(p->oset) );
DSET_delete( p->oset );
RETURN(-1);
}
tross_Copy_History( p->gpar, p->oset );
tross_Make_History( PROG_NAME, argc, argv, p->oset );
if ( DSET_write( p->oset ) != True )
{
fprintf( stderr, "** failed to write dataset '%s', exiting...\n",
opts->oset_file );
RETURN(-1);
}
RETURN(0);
}
/*----------------------------------------------------------------------
* s2v_nodes2volume - create an AFNI dataset from the node_list
* and map type (subject to any mask).
*
* inputs:
* N - xyz coordinate node list for surfaces
* p - param_t struct
* map - the surface to dataset mapping type
*
* return AFNI dataset - on success
* NULL - on failure
*----------------------------------------------------------------------
*/
THD_3dim_dataset * s2v_nodes2volume( node_list_t * N, param_t * p,
s2v_opts_t * sopt )
{
THD_3dim_dataset * dout;
THD_fvec3 * pary;
double * ddata;
void * vdata = NULL;
int * cdata;
float fac;
int nvox, dsize, valid;
int sub;
ENTRY("s2v_nodes2volume");
if ( N == NULL || ! ISVALID_DSET(p->gpar) )
{
fprintf(stderr,"** s2v_nodes2volume: bad params (%p,%p)\n", N, p->gpar);
RETURN(NULL);
}
dout = EDIT_empty_copy( p->gpar );
if ( ! ISVALID_3DIM_DATASET( dout ) )
{
fprintf( stderr, "** failed EDIT_empty_copy()\n" );
RETURN(NULL);
}
/* output dataset will have nsubs sub-bricks */
EDIT_dset_items( dout, ADN_nvals, p->nsubs,
ADN_type, HEAD_FUNC_TYPE,
ADN_func_type, FUNC_BUCK_TYPE,
ADN_datum_all, sopt->datum,
ADN_ntt, 0,
ADN_none );
if ( sopt->debug > 1 )
fprintf( stderr, "++ creating dataset '%s' of type '%s', nvals = %d\n",
DSET_HEADNAME(dout), MRI_TYPE_name[sopt->datum], p->nsubs );
nvox = DSET_NVOX(p->gpar);
/* allocate a computational array of doubles */
if ( (ddata = (double *)malloc(nvox * sizeof(double))) == NULL )
{
fprintf( stderr, "** n2v: failed to allocate %d doubles\n", nvox );
DSET_delete( dout );
RETURN(NULL);
}
/* allocate an array of ints for functional counting */
if ( (cdata = (int *)malloc(nvox * sizeof(int))) == NULL )
{
fprintf( stderr, "** n2v: failed to allocate %d ints\n", nvox );
DSET_delete( dout ); free( ddata );
RETURN(NULL);
}
/* allocate space for the point list */
if ( (pary = (THD_fvec3 *)malloc(sopt->f_steps*sizeof(THD_fvec3))) == NULL )
{
fprintf(stderr,"** n2v: cannot allocate %d THD_fvec3s\n",sopt->f_steps);
DSET_delete( dout ); free( ddata ); free( cdata );
RETURN(NULL);
}
dsize = mri_datum_size(sopt->datum);
/* create the sub-brick data for output */
for ( sub = 0; sub < p->nsubs; sub++ )
{
valid = 0;
/* need new data after each sub-brick is attached to the dset */
if ( (vdata = malloc( nvox * dsize )) == NULL )
{
fprintf( stderr, "** failed to allocate %d bytes for vdata\n",
nvox * dsize );
DSET_delete( dout ); free( ddata ); free( cdata );
RETURN(NULL);
}
/* init data to zeros */
memset(cdata, 0, nvox*sizeof(int));
memset(ddata, 0, nvox*sizeof(double));
/* set node_list data */
if ( set_node_list_data( N, p, sopt, sub ) != 0 )
{
DSET_delete(dout);
free(ddata);
free(cdata);
free(vdata);
RETURN(NULL);
}
if ( compute_results( p, N, sopt, ddata, cdata, pary ) == 0 )
valid = 1;
if ( ! valid ) /* then clean up memory */
{
free( ddata );
free( vdata );
DSET_delete( dout );
RETURN(NULL);
}
/* convert to output data type */
fac = 0.0;
if ( sopt->debug > 1 )
fprintf(stderr,"++ sub-brick %d, integral_doubles = %d\n",
sub, integral_doubles( ddata, nvox ));
/* for int output and scaling, fac is maxval/maxabs */
if ( MRI_IS_INT_TYPE( sopt->datum ) && !sopt->noscale )
{
float amax = MCW_vol_amax( nvox, 1, 1, MRI_double, ddata );
if ( amax > MRI_TYPE_maxval[sopt->datum] ) /* we must scale */
fac = MRI_TYPE_maxval[sopt->datum]/amax;
else if ( integral_doubles( ddata, nvox ) ) /* don't need scale */
fac = 0.0;
else if ( amax != 0.0 )
fac = MRI_TYPE_maxval[sopt->datum]/amax;
if ( sopt->debug > 1 )
fprintf(stderr,"++ fac = %f, amax = %f \n", fac, amax);
}
EDIT_coerce_scale_type(nvox, fac, MRI_double,ddata, sopt->datum,vdata);
/* attach the final vdata to the dataset */
EDIT_substitute_brick( dout, sub, sopt->datum, vdata );
DSET_BRICK_FACTOR( dout, sub ) = (fac != 0.0) ? 1.0/fac : 0.0;
}
if (sopt->debug > 0)
fprintf(stderr,"++ %d sub-brick(s) computed\n", p->nsubs);
free(ddata);
free(cdata);
free(pary);
RETURN(dout);
}
/*----------------------------------------------------------------------
* compute_results - fill ddata with results, cdata with count
*
* return 0 : success
* else : failure
*----------------------------------------------------------------------
*/
int compute_results( param_t * p, node_list_t * N, s2v_opts_t * sopt,
double * ddata, int * cdata, THD_fvec3 * pary )
{
THD_fvec3 p1, pn;
float dist, min_dist, max_dist;
int nindex, node;
int oobc;
ENTRY("compute_results");
if ( !p || !N || !sopt || !ddata || !cdata )
{
fprintf(stderr,"** cr: bad params (%p,%p,%p,%p,%p)\n",
p, N, sopt, ddata, cdata);
RETURN(-1);
}
min_dist = 9999.9;
max_dist = -1.0;
oobc = 0; /* init out-of-bounds counter */
for ( nindex = 0; nindex < N->ilen; nindex++ )
{
node = N->ilist[nindex];
if ( node < 0 || node >= N->nnodes )
{
fprintf(stderr,"** node %d (%d) out-of-range [0,%d]\n",
nindex, node, N->nnodes);
RETURN(-1);
}
/* initiate the endpoints */
p1 = N->nodes[node];
if ( N->depth > 1 )
pn = N->nodes[node+N->nnodes];
else
pn = p1;
if ( ! sopt->sxyz_ori_gpar ) /* if orient is not as gpar it is Dicom */
{
p1 = THD_dicomm_to_3dmm(p->gpar, p1);
pn = THD_dicomm_to_3dmm(p->gpar, pn);
}
if ( sopt->debug > 0 && sopt->dnode == node )
{
fprintf(stderr,"-- debug node: %d\n", node );
fprintf(stderr,"-- orig endpts: (%f, %f, %f)\n"
" (%f, %f, %f)\n",
p1.xyz[0], p1.xyz[1], p1.xyz[2],
pn.xyz[0], pn.xyz[1], pn.xyz[2]);
}
adjust_endpts( sopt, &p1, &pn ); /* per user options */
/* if both points are outside gpar box, skip the node */
if ( f3mm_out_of_bounds( &p1, &p->f3mm_min, &p->f3mm_max ) &&
f3mm_out_of_bounds( &pn, &p->f3mm_min, &p->f3mm_max ) )
{
oobc++;
continue;
}
dist = dist_f3mm( &p1, &pn );
if ( dist < min_dist ) min_dist = dist;
if ( dist > max_dist ) max_dist = dist;
make_point_list( pary, &p1, &pn, sopt );
/* do all the work to insert data */
if ( insert_list(N, p, sopt, pary, nindex, ddata, cdata) )
RETURN(-1);
}
if ( sopt->debug > 0 )
fprintf(stderr, "-- (min_dist, max_dist) = (%f, %f)\n"
" %d of %d nodes were out of bounds\n",
min_dist, max_dist, oobc, N->ilen);
if ( final_computations(ddata, cdata, sopt, DSET_NVOX(p->gpar)) )
RETURN(-1);
RETURN(0);
}
/*----------------------------------------------------------------------
* final_computations - perform any last computation over the dataset
*----------------------------------------------------------------------
*/
int final_computations(double *ddata, int *cdata, s2v_opts_t *sopt, int nvox)
{
int index;
ENTRY("final_computations");
if ( (sopt->debug > 1) && (sopt->dvox >= 0) && (sopt->dvox < nvox) )
fprintf(stderr,"++ final: voxel %d, from values (%d, %f) ",
sopt->dvox,cdata[sopt->dvox],ddata[sopt->dvox]);
switch ( sopt->map )
{
default:
fprintf(stderr,"** fc: mapping %d not ready\n", sopt->map );
RETURN(-1);
case E_SMAP_AVE:
/* we have sum, divide for average */
for ( index = 0; index < nvox; index++ )
if ( cdata[index] > 0 )
ddata[index] /= cdata[index];
break;
case E_SMAP_COUNT:
case E_SMAP_MAX_ABS:
case E_SMAP_MAX:
case E_SMAP_MIN:
case E_SMAP_MASK:
case E_SMAP_MASK2:
break;
}
if ( (sopt->debug > 1) && (sopt->dvox >= 0) && (sopt->dvox < nvox) )
fprintf(stderr,"to values (%d, %f)\n",
cdata[sopt->dvox],ddata[sopt->dvox]);
RETURN(0);
}
/*----------------------------------------------------------------------
* insert_list - given xyz list, add to ddata/cdata
*
* for each point in list
* get voxel index
* if out of mask, skip
* if undesired voxel repeat, skip
* apply mapping function to insert data
*
* return 0 on success
* < 0 on error
*----------------------------------------------------------------------
*/
int insert_list( node_list_t * N, param_t * p, s2v_opts_t * sopt,
THD_fvec3 * pary, int nindex, double * ddata, int * cdata )
{
THD_ivec3 i3;
int vindex, prev_vind;
int step, node;
int nx, ny, debug;
ENTRY("insert_list");
if ( !N || !p || !sopt || !pary || nindex < 0 || !ddata || !cdata )
{
fprintf(stderr,"** ID params (%p,%p,%p,%p,%d,%p,%p)\n",
N, p, sopt, pary, nindex, ddata, cdata);
RETURN(-1);
}
node = N->ilist[nindex];
debug = sopt->debug && node == sopt->dnode;
if ( debug ) /* show xyz coords? */
fprintf(stderr,"++ point list for N[%d] = %d:\n", nindex, node);
nx = DSET_NX(p->gpar);
ny = DSET_NY(p->gpar);
prev_vind = -1;
for ( step = 0; step < sopt->f_steps; step++ )
{
if ( f3mm_out_of_bounds( pary+step, &p->f3mm_min, &p->f3mm_max ) )
continue;
i3 = THD_3dmm_to_3dind(p->gpar, pary[step]);
vindex = i3.ijk[0] + nx * (i3.ijk[1] + ny * i3.ijk[2]);
if ( debug )
fprintf(stderr," %3d, vox %d, [%3d,%3d,%3d]: %f, %f, %f",
step, vindex, i3.ijk[0], i3.ijk[1], i3.ijk[2],
pary[step].xyz[0], pary[step].xyz[1], pary[step].xyz[2]);
if ( sopt->cmask && !sopt->cmask[vindex] )
{
if ( debug ) fprintf(stderr, " : skip (mask)\n");
continue;
}
if ( (vindex == prev_vind) && (sopt->f_index == S2V_F_INDEX_VOXEL) )
{
if ( debug ) fprintf(stderr, " : skip (repeat voxel)\n");
continue;
}
/* hey, we finally have a new voxel to add */
prev_vind = vindex;
if ( debug )
fprintf( stderr, " : (old) %d %f", cdata[vindex],ddata[vindex]);
if (insert_value(sopt, ddata, cdata, vindex, node, N->fdata[nindex]))
RETURN(-1);
if ( debug )
fprintf( stderr, " : (new) %d %f\n", cdata[vindex],ddata[vindex]);
}
RETURN(0);
}
/*----------------------------------------------------------------------
* insert_value - insert value based on user function
*
* return 0 on success
* < 0 on error
*----------------------------------------------------------------------
*/
int insert_value(s2v_opts_t * sopt, double *dv, int *cv, int vox, int node,
float value)
{
ENTRY("insert_value");
if ( !sopt || !dv || !cv || vox < 0 )
{
fprintf(stderr, "** IV: bad params (%p,%p,%p,%d)\n",sopt,dv,cv,vox);
RETURN(-1);
}
if ( (sopt->debug > 1) && (vox == sopt->dvox) )
fprintf(stderr,"++ voxel %d, node %d, values (%d, %f) ",
vox, node, cv[vox], dv[vox]);
switch ( sopt->map )
{
default:
fprintf(stderr,"** IV: mapping %d not ready\n", sopt->map );
RETURN(-1);
case E_SMAP_AVE:
if ( cv[vox] == 0 )
dv[vox] = value;
else
dv[vox] += value; /* divide by count later */
break;
case E_SMAP_COUNT:
if ( cv[vox] == 0 )
dv[vox] = 1.0;
else
dv[vox]++;
break;
case E_SMAP_MAX:
if ( cv[vox] == 0 )
dv[vox] = value;
else if (value > dv[vox])
dv[vox] = value;
break;
case E_SMAP_MAX_ABS:
if ( cv[vox] == 0 )
dv[vox] = value;
else if (fabs(value) > fabs(dv[vox]))
dv[vox] = value;
break;
case E_SMAP_MIN:
if ( cv[vox] == 0 )
dv[vox] = value;
else if (value < dv[vox])
dv[vox] = value;
break;
case E_SMAP_MASK:
case E_SMAP_MASK2:
dv[vox] = 1.0;
break;
}
cv[vox]++; /* in any case, increment the counter */
if ( (sopt->debug > 1) && (vox == sopt->dvox) )
fprintf(stderr,"to (%d, %f)\n",cv[vox],dv[vox]);
RETURN(0);
}
/*----------------------------------------------------------------------
* make_point_list - from endpoints and steps
*
* return 0 on success
* < 0 on error
*----------------------------------------------------------------------
*/
int make_point_list( THD_fvec3 * list, THD_fvec3 * p1, THD_fvec3 * pn,
s2v_opts_t * sopt )
{
float rat1, ratn;
int step;
ENTRY("make_point_list");
if ( !list || !p1 || !pn || !sopt )
{
fprintf(stderr,"** mpl: bad params (%p,%p,%p,%p)\n", list,p1,pn,sopt);
RETURN(-1);
}
list[0] = *p1;
for ( step = 1; step < sopt->f_steps; step++ )
{
ratn = step / (sopt->f_steps - 1.0);
rat1 = 1.0 - ratn;
list[step].xyz[0] = rat1 * p1->xyz[0] + ratn * pn->xyz[0];
list[step].xyz[1] = rat1 * p1->xyz[1] + ratn * pn->xyz[1];
list[step].xyz[2] = rat1 * p1->xyz[2] + ratn * pn->xyz[2];
}
RETURN(0);
}
/*----------------------------------------------------------------------
* adjust_endpts - adjust endpoints for map and options
*
* return 0 on success
* < 0 on error
*----------------------------------------------------------------------
*/
int adjust_endpts( s2v_opts_t * sopt, THD_fvec3 * p1, THD_fvec3 * pn )
{
THD_fvec3 f3_diff;
float dist, factor;
ENTRY("adjust_endpts");
if ( !sopt || !p1 || !pn )
{
fprintf(stderr,"** ae: invalid params (%p,%p,%p)\n", sopt, p1, pn);
RETURN(-1);
}
/* first, get the difference, and distance */
f3_diff.xyz[0] = pn->xyz[0] - p1->xyz[0];
f3_diff.xyz[1] = pn->xyz[1] - p1->xyz[1];
f3_diff.xyz[2] = pn->xyz[2] - p1->xyz[2];
dist = dist_f3mm( p1, pn );
if ( (sopt->f_p1_fr != 0.0) || (sopt->f_p1_mm != 0.0) )
{
if ( sopt->f_p1_fr != 0.0 ) /* what the heck, choose fr if both */
factor = sopt->f_p1_fr;
else
factor = (dist == 0.0) ? 0.0 : sopt->f_p1_mm / dist;
p1->xyz[0] += factor * f3_diff.xyz[0];
p1->xyz[1] += factor * f3_diff.xyz[1];
p1->xyz[2] += factor * f3_diff.xyz[2];
}
if ( (sopt->f_pn_fr != 0.0) || (sopt->f_pn_mm != 0.0) )
{
if ( sopt->f_pn_fr != 0.0 )
factor = sopt->f_pn_fr;
else
factor = (dist == 0.0) ? 0.0 : sopt->f_pn_mm / dist;
pn->xyz[0] += factor * f3_diff.xyz[0];
pn->xyz[1] += factor * f3_diff.xyz[1];
pn->xyz[2] += factor * f3_diff.xyz[2];
}
switch ( sopt->map )
{
default:
fprintf(stderr,"** ae: mapping %d not ready\n", sopt->map );
RETURN(-1);
case E_SMAP_AVE:
case E_SMAP_COUNT:
case E_SMAP_MAX_ABS:
case E_SMAP_MAX:
case E_SMAP_MIN:
case E_SMAP_MASK:
case E_SMAP_MASK2:
break;
}
RETURN(0);
}
/*----------------------------------------------------------------------
* set_node_list_data - fill fdata in node list
*
* For the given column, fill N->ilist with data from p->sdata_im.
*----------------------------------------------------------------------
*/
int set_node_list_data ( node_list_t *N, param_t *p, s2v_opts_t *sopt, int col )
{
float * fp, fval;
int c, lposn;
ENTRY("set_node_list_data");
if ( !N || !p || !sopt || col < 0 )
{
fprintf(stderr,"** snld: bad params (%p,%p,%p,%d)\n", N, p, sopt, col);
RETURN(-1);
}
if ( sopt->debug > 1 )
fprintf(stderr, "-- setting fdata for column %d\n", col);
/* if sdata_im does not exist, fill as mask and return */
if ( !p->sdata_im )
{
fval = 1.0; /* init to mask value */
/* if we have a PARSER_code (presumably without symbols), use it */
if ( p->parser.pcode )
{
for ( c = 0; c < 26; c++ )
p->parser.atoz[c] = 0.0;
fval = PARSER_evaluate_one( p->parser.pcode, p->parser.atoz );
}
for ( c = 0; c < N->ilen; c++ )
N->fdata[c] = fval;
RETURN(0);
}
/* else sdata exists: check column, and copy data from sdata_im */
if ( col > (p->sdata_im->nx - 2) && !p->parser.pcode )
{
fprintf(stderr,"** snld error: col > nx-2 (%d > %d)\n",
col, p->sdata_im->nx-2);
RETURN(-1);
}
else if ( p->sdata_im->ny < N->ilen )
{
fprintf(stderr,"** snld error: ny < ilen (%d < %d)\n",
p->sdata_im->ny, N->ilen);
RETURN(-1);
}
else if ( !N->fdata )
{
fprintf(stderr,"** snld error: missing idata\n");
RETURN(-1);
}
else if ( p->parser.pcode && col != 0 ) /* let's be safe */
{
fprintf(stderr,"** snld error: cannot use parser with col = %d\n", col);
RETURN(-1);
}
/* hmmmm, we're still missing something... oh yes, data! */
fp = MRI_FLOAT_PTR( p->sdata_im ) + col+1; /* offset by column number */
for ( c = 0; c < N->ilen; c++ )
{
if ( p->parser.pcode )
{
/* fill atoz with surface node data */
for ( lposn = 0; lposn < p->parser.max_sym; lposn++ )
p->parser.atoz[lposn] = fp[lposn];
N->fdata[c] = PARSER_evaluate_one(p->parser.pcode, p->parser.atoz);
}
else
N->fdata[c] = *fp;
fp += p->sdata_im->nx;
}
RETURN(0);
}
/*----------------------------------------------------------------------
* init_node_list - from surfaces or from sxyz file
*
* Fill nodes with coordinate data.
*----------------------------------------------------------------------
*/
int init_node_list (opts_t *opts, param_t *p, s2v_opts_t *sopt, node_list_t *N)
{
int nsurf, rv;
ENTRY("init_node_list");
if ( opts == NULL || p == NULL || sopt == NULL || N == NULL )
{
fprintf(stderr, "** inl - bad params (%p,%p,%p,%p)\n",opts,p,sopt,N);
RETURN(-1);
}
memset(N,0,sizeof(*N)); /* first, clear it out - have care w/pointers */
N->nodes = NULL; N->fdata = NULL; N->ilist = NULL;
nsurf = SUMAg_N_DOv; /* verify, as long as there is no sxyz_im */
if ( !p->sxyz_im && ((nsurf < 1) || (nsurf > S2V_MAX_SURFS)) )
{
fprintf(stderr,"** inl: SUMAg_N_DOv has invalid value of %d\n", nsurf);
RETURN(-1);
}
#if 0 /* remove requrement of 2 surfaces for functions [v3.0] */
if ( sopt->map == E_SMAP_MASK )
nsurf = 1;
else
{
/* do a quick check before slow reading action */
if ( ! p->sxyz_im && (SUMAg_N_DOv < 2) )
{
fprintf(stderr,"** function '%s' requires 2 surfaces\n",
gs2v_map_names[sopt->map]);
RETURN(-1);
}
nsurf = 2;
}
#endif
if ( p->sxyz_im )
rv = sxyz_1D_to_nlist( opts, sopt, p, N, &nsurf );
else
rv = surf_to_node_list( sopt, N, nsurf );
if ( sopt->debug > 1 && rv == 0 ) /* errors reported separately */
fprintf(stderr,"++ node list initialized\n");
RETURN(rv);
}
/*----------------------------------------------------------------------
* sxyz_1D_to_nlist - nodes from sxyz 1D to node_list
*
*----------------------------------------------------------------------
*/
int sxyz_1D_to_nlist(opts_t * opts, s2v_opts_t * sopt, param_t * p,
node_list_t * N, int * nsurf)
{
THD_fvec3 * fvp;
float * fp;
int c, sc;
ENTRY("sxyz_1D_to_nlist");
if ( !sopt || !p || !N )
{
fprintf(stderr,"** sxyz2nl: bad params (%p,%p,%p)\n",sopt,p,N);
RETURN(-1);
}
if ( !p->sxyz_im )
{
fprintf(stderr,"** missing sxyz_im for sxyz surf '%s'\n",
opts->surf_xyz_1D_file);
RETURN(-1);
}
*nsurf = p->sxyz_im->nx / 3;
if ( p->sxyz_im->nx != 3 * *nsurf )
{
fprintf(stderr,"** sxyz surf '%s' has %d columns (%d expected)\n",
opts->surf_xyz_1D_file, p->sxyz_im->nx, 3**nsurf);
RETURN(-1);
}
else if ( p->sxyz_im->ny <= 0 )
{
fprintf(stderr,"** sxyz surf '%s': bad sxyz dimensions (%d,%d)\n",
opts->surf_xyz_1D_file, p->sxyz_im->nx, p->sxyz_im->ny);
RETURN(-1);
}
N->depth = *nsurf;
N->nnodes = p->sxyz_im->ny;
N->nodes = (THD_fvec3 *)malloc(N->depth*N->nnodes*sizeof(THD_fvec3));
if ( N->nodes == NULL )
{
fprintf(stderr,"** failed to allocate %dx%d THD_fvec3's for nodes\n",
N->nnodes, N->depth);
RETURN(-1);
}
fvp = N->nodes;
for ( sc = 0; sc < *nsurf; sc++ )
{
fp = MRI_FLOAT_PTR( p->sxyz_im ) + sc * 3; /* offset for surf */
for ( c = 0; c < N->nnodes; c++ )
{
fvp->xyz[0] = fp[0];
fvp->xyz[1] = fp[1];
fvp->xyz[2] = fp[2];
fp += p->sxyz_im->nx;
fvp++;
}
}
if ( sopt->debug > 0 )
{
fprintf(stderr, "++ sxyz_1D nodes from '%s': nxyz = %d, nsurf = %d\n",
opts->surf_xyz_1D_file, N->nnodes, N->depth);
if ( sopt->dnode >= 0 && sopt->dnode <= p->sxyz_im->ny )
{
fprintf(stderr," debug node (%d) loc", sopt->dnode);
for ( sc = 0; sc < *nsurf; sc++ )
fprintf(stderr," : (%f, %f, %f)",
N->nodes[sopt->dnode + sc*N->nnodes].xyz[0],
N->nodes[sopt->dnode + sc*N->nnodes].xyz[1],
N->nodes[sopt->dnode + sc*N->nnodes].xyz[2]);
fputc('\n', stderr);
}
}
RETURN(0);
}
/*----------------------------------------------------------------------
* surf_to_node_list - create node list for mask mapping
*
* Allocate space for node list(s), and fill with xyz coordinates.
*
* return -1 : on error
* 0 : on success
*----------------------------------------------------------------------
*/
int surf_to_node_list ( s2v_opts_t * sopt, node_list_t * N, int nsurf )
{
SUMA_SurfaceObject ** so;
THD_fvec3 * fvp;
float * fp;
int rv, nindex, sindex;
ENTRY("surf_to_node_list");
if ( sopt == NULL || N == NULL || nsurf < 0 || nsurf > 2 )
{
fprintf( stderr, "** anl: bad params (%p,%p,%d)\n",
sopt, N, nsurf );
RETURN(-1);
}
/* create a temporary list of surface pointers */
so = (SUMA_SurfaceObject **)calloc(nsurf, sizeof(SUMA_SurfaceObject *));
if ( so == NULL )
{
fprintf( stderr, "** anl: failed to alloc %d surf pointers\n", nsurf );
RETURN(-1);
}
if ( (rv = get_mappable_surfs( so, nsurf, sopt->debug )) != nsurf )
{
fprintf( stderr, "** error: found %d (of %d) mappable surfaces\n",
rv, nsurf );
RETURN(-1);
}
/* fill node list struct */
N->depth = nsurf;
N->nnodes = so[0]->N_Node;
N->nodes = (THD_fvec3 *)malloc(N->depth * N->nnodes * sizeof(THD_fvec3));
if ( N->nodes == NULL )
{
fprintf( stderr, "** cnlm: failed to allocate %d THD_fvec3 structs\n",
N->depth * N->nnodes );
free(so);
RETURN(-1);
}
/* copy the xyz coordinates for each node */
fvp = N->nodes; /* linear coverage of all nodes */
for ( sindex = 0; sindex < N->depth; sindex++ )
{
if ( so[sindex]->N_Node != N->nnodes )
{
fprintf(stderr, "** surf #%d (%s) has %d nodes (but expected %d)\n",
sindex,
so[sindex]->Label ? so[sindex]->Label : "<unnamed>",
so[sindex]->N_Node, N->nnodes );
free( N->nodes ); N->nodes = NULL;
free(so);
RETURN(-1);
}
for ( nindex = 0, fp = so[sindex]->NodeList;
nindex < N->nnodes;
nindex++, fp += 3 )
{
memcpy( fvp->xyz, fp, 3*sizeof(float) );
fvp++;
}
}
if ( sopt->debug > 1 )
fprintf( stderr, "++ allocated %d x %d (x %d) node list\n",
N->depth, N->nnodes, (int)sizeof(THD_fvec3) );
free(so);
RETURN(0);
}
/*----------------------------------------------------------------------
* get_mappable_surfs - return mappable surface objects
*
* return the number of surfaces found
*----------------------------------------------------------------------
*/
int get_mappable_surfs( SUMA_SurfaceObject ** slist, int how_many, int debug )
{
SUMA_SurfaceObject * so;
int count, socount = 0;
ENTRY("get_mappable_surfts");
if ( slist == NULL )
{
fprintf( stderr, "** gms: missing slist!\n" );
RETURN(-1);
}
for ( count = 0; count < SUMAg_N_DOv; count++ )
{
if ( ! SUMA_isSO(SUMAg_DOv[count]) )
continue;
so = (SUMA_SurfaceObject *)SUMAg_DOv[count].OP;
if ( ! so->AnatCorrect )
{
if ( debug )
fprintf(stderr,"-- surf #%d '%s', anat not correct, skipping\n",
socount, CHECK_NULL_STR(so->Label));
if ( debug > 1 )
fprintf(stderr,"** consider adding the following to the "
"surface definition in the spec file:\n"
" Anatomical = Y\n");
continue;
}
if ( debug > 1 )
{
fprintf( stderr, "\n---------- surface #%d '%s' -----------\n",
socount, CHECK_NULL_STR(so->Label) );
SUMA_Print_Surface_Object( so, stderr );
}
if ( socount < how_many ) /* store a good surface */
slist[socount] = so;
socount++;
}
if ( debug > 1 )
fprintf( stderr, "++ found %d mappable surfaces\n", socount );
RETURN(socount);
}
/*----------------------------------------------------------------------
* set_smap_opts - given options and mapping function, fill struct
*
* return 0 : success
* -1 : error condition
*----------------------------------------------------------------------
*/
int set_smap_opts( opts_t * opts, param_t * p, s2v_opts_t * sopt )
{
int nsurf;
ENTRY("set_smap_opts");
memset( sopt, 0, sizeof(*sopt) ); /* clear the sopt struct */
sopt->cmask = NULL;
if ( (sopt->map = check_map_func( opts->map_str )) == E_SMAP_INVALID )
RETURN(-1);
sopt->datum = check_datum_type(opts->datum_str, DSET_BRICK_TYPE(p->gpar,0));
if (sopt->datum < 0)
RETURN(-1);
if ( opts->noscale == 1 )
sopt->noscale = 1;
sopt->debug = opts->debug; /* for output in library functions */
sopt->dnode = opts->dnode;
sopt->dvox = opts->dvox;
sopt->sxyz_ori_gpar = opts->sxyz_ori_gpar;
sopt->cmask = p->cmask;
if ( opts->f_steps < S2V_F_STEPS_MIN ) /* minimum is 1 */
{
/* if f_steps was not given, init to the number of surfaces - v3.6 */
for ( nsurf = 0; nsurf < S2V_MAX_SURFS && opts->snames[nsurf]; nsurf++ )
;
if ( nsurf <= 0 )
{
fprintf(stderr,"** error: sso: no input surfaces\n");
RETURN(-1);
}
sopt->f_steps = nsurf;
}
else
sopt->f_steps = opts->f_steps;
sopt->f_index = S2V_F_INDEX_VOXEL;
if ( (opts->f_index_str != NULL) &&
( !strncmp(opts->f_index_str, "point", 5) ||
!strncmp(opts->f_index_str, "node", 4) ) ) /* preserve old use */
sopt->f_index = S2V_F_INDEX_POINT;
sopt->f_p1_fr = opts->f_p1_fr; /* copy fractions & distances */
sopt->f_pn_fr = opts->f_pn_fr;
sopt->f_p1_mm = opts->f_p1_mm;
sopt->f_pn_mm = opts->f_pn_mm;
switch (sopt->map)
{
default:
case E_SMAP_AVE:
case E_SMAP_MAX:
case E_SMAP_MIN:
case E_SMAP_MAX_ABS:
break;
case E_SMAP_COUNT:
case E_SMAP_MASK:
case E_SMAP_MASK2:
sopt->noscale = 1;
break;
}
if ( opts->debug > 1 )
disp_s2v_opts_t( "++ s2v_opts_set :", sopt );
RETURN(0);
}
/*----------------------------------------------------------------------
* free memory, close output file
*----------------------------------------------------------------------
*/
int final_clean_up ( node_list_t * N )
{
ENTRY("final_clean_up");
if ( ( SUMAg_DOv != NULL ) &&
( SUMA_Free_Displayable_Object_Vect(SUMAg_DOv, SUMAg_N_DOv) == 0 ) )
fprintf(stderr, "** failed SUMA_Free_Displayable_Object_Vect()\n" );
if ( ( SUMAg_SVv != NULL ) &&
( SUMA_Free_SurfaceViewer_Struct_Vect(SUMAg_SVv, SUMAg_N_SVv) == 0 ) )
fprintf( stderr, "** failed SUMA_Free_SurfaceViewer_Struct_Vect()\n" );
if ( ( SUMAg_CF != NULL ) && ( SUMA_Free_CommonFields(SUMAg_CF) == 0 ) )
fprintf( stderr, "** failed SUMA_Free_CommonFields()\n" );
if ( N )
{
if ( N->nodes ) free( N->nodes );
if ( N->fdata ) free( N->fdata );
if ( N->ilist ) free( N->ilist );
}
RETURN(0);
}
/*----------------------------------------------------------------------
*
*----------------------------------------------------------------------
*/
int read_surf_files ( opts_t * opts, param_t * p, SUMA_SurfSpecFile * spec,
s2v_opts_t * sopt, node_list_t * N )
{
int rv;
ENTRY("read_surf_files");
/* get surface coordinates */
if ( opts->spec_file )
{
if ( (rv = fill_SUMA_structs(opts, spec)) != 0 )
RETURN(rv);
}
else if ( opts->surf_xyz_1D_file )
{
if ( (rv = read_sxyz_1D( opts, p )) != 0 )
RETURN(rv);
}
else
{
fprintf(stderr,"** missing spec or sdata file, exiting...\n");
RETURN(-1);
}
/* assign node list nodes (xyz coords) */
if ( (rv = init_node_list(opts, p, sopt, N)) != 0 )
RETURN(rv);
/* check status and allocate memory in node list */
if ( (rv = fill_node_list(opts, p, N)) != 0 )
RETURN(rv);
if ( (rv = verify_parser_expr(opts, p)) != 0 )
RETURN(rv);
RETURN(0);
}
/*----------------------------------------------------------------------
* read_sxyz_file - read surf_xyz_1D_file
*----------------------------------------------------------------------
*/
int read_sxyz_1D ( opts_t * opts, param_t * p )
{
MRI_IMAGE * im;
ENTRY("read_sxyz_1D");
if ( !opts || !p )
{
fprintf(stderr,"** rsxyz1D: bad params (%p,%p)\n", opts, p);
RETURN(-1);
}
if ( (im = mri_read_1D(opts->surf_xyz_1D_file)) == NULL )
{
fprintf(stderr,"** failed to read sxyz 1D file '%s'\n",
opts->surf_xyz_1D_file);
RETURN(-1);
}
if ( im->nx < 1 || im->ny < 3 ) /* must have xyz coords, at least */
{
fprintf(stderr,"** bad sxyz file '%s'?\n", opts->surf_xyz_1D_file);
RETURN(-1);
}
/* transpose list to node-major order, and lose temp image */
p->sxyz_im = mri_transpose(im);
mri_free(im);
if ( opts->debug > 0 )
fprintf(stderr,"++ read 1D xyz surface file '%s' (nx = %d, ny = %d)\n",
opts->surf_xyz_1D_file, p->sxyz_im->nx, p->sxyz_im->ny );
RETURN(0);
}
/*----------------------------------------------------------------------
* fill the node_list data fields
*
* ilist gets index list, if any
* fdata gets allocation for data
* ilen is the length of the two lists
*----------------------------------------------------------------------
*/
int fill_node_list ( opts_t * opts, param_t * p, node_list_t * N )
{
int rv;
ENTRY("fill_node_list");
if ( !opts || !p || !N )
{
fprintf(stderr,"** fill_node_list: bad params (%p,%p,%p)\n",opts,p,N);
RETURN(-1);
}
p->nsubs = 1; /* unless we get additional surface data */
if ( opts->sdata_file_1D )
{
if ( (rv = sdata_from_1D( opts, p, N )) != 0 )
RETURN(rv);
if ( ! p->parser.pcode )
p->nsubs = p->sdata_im->nx - 1;
}
else
{
if ( (rv = sdata_from_default( N )) != 0 )
RETURN(rv);
}
if ( (rv = verify_node_list( N )) != 0 )
RETURN(rv);
if ( opts->debug > 1 )
disp_node_list_t( "++ node list filled: ", N );
RETURN(0);
}
/*----------------------------------------------------------------------
* verify that the parser expression has sufficient surface data values
*----------------------------------------------------------------------
*/
int verify_parser_expr ( opts_t * opts, param_t * p )
{
int max_used;
ENTRY("verify_parser_expr");
if ( !opts || !p )
{
fprintf(stderr,"** vpe: invalid params (%p,%p)\n", opts, p);
RETURN(-1);
}
/* if no parser code, there is nothing to do */
if ( ! p->parser.pcode )
RETURN(0);
for ( max_used = 25; max_used >= 0; max_used-- )
if ( p->parser.has_sym[max_used] )
break;
max_used++; /* this is the number of surface values needed */
p->parser.max_sym = max_used;
/* if the expression is not constant, we need some data */
if ( max_used > 0 )
{
if ( !p->sdata_im )
{
fprintf(stderr, "** parser expression requires surface data\n"
" (see '-sdata_1D')\n");
RETURN(-1);
}
else if ( max_used > p->sdata_im->nx - 1 )
{
fprintf(stderr,
"** error: not enough surface values for expression\n"
" svals = %d, exp_vals = %d, expr = '%s'\n",
p->sdata_im->nx - 1, max_used, opts->data_expr);
RETURN(-1);
}
}
if ( opts->debug > 1 )
fprintf(stderr,"-- surf_vals = %d, expr_vals = %d\n",
p->sdata_im ? (p->sdata_im->nx - 1) : 0, max_used);
RETURN(0);
}
/*----------------------------------------------------------------------
* verify that the node list makes sense
*----------------------------------------------------------------------
*/
int verify_node_list ( node_list_t * N )
{
int icount, errs = 0;
ENTRY("verify_node_list");
if ( !N )
{
fprintf(stderr, "** vnl - no node list\n" );
RETURN(-1);
}
if ( !N->nodes || !N->ilist )
{
fprintf(stderr,"** missing nodes or ilist\n" );
errs++;
}
if ( N->depth < 1 || N->nnodes < 1 || N->ilen < 1 )
{
fprintf(stderr,"** invalid depth, nnodes or ilen" );
errs++;
}
if ( errs )
{
disp_node_list_t("** invalid data : ", N );
RETURN(-errs);
}
/* now check that the indices are within nnodes range */
for ( icount = 0; icount < N->ilen; icount++ )
if ( N->ilist[icount] < 0 || N->ilist[icount] >= N->nnodes )
{
fprintf(stderr,"** surf data index number %d is out of range:\n"
" index = %d, range is [%d,%d]\n",
icount, N->ilist[icount], 0, N->nnodes-1);
RETURN(-10);
}
/* node_list is okay, so make space for the actual data */
if ( (N->fdata = (float *)malloc(N->ilen*sizeof(float))) == NULL )
{
fprintf(stderr,"** vnl: failed to allocate %d floats\n",N->ilen);
free(N->ilist);
RETURN(-20);
}
RETURN(0);
}
/*----------------------------------------------------------------------
* fill ilist with default indices (i[k] = k)
*----------------------------------------------------------------------
*/
int sdata_from_default ( node_list_t * N )
{
int c;
ENTRY("sdata_from_default");
if ( !N )
RETURN(-1);
if ( (N->ilist = (int *)malloc(N->nnodes * sizeof(int))) == NULL )
{
fprintf(stderr,"** failed to allocate %d ints for ilist\n",N->nnodes);
RETURN(-1);
}
N->ilen = N->nnodes;
/* now fill ilist with default node indices */
for ( c = 0; c < N->ilen; c++ )
N->ilist[c] = c;
RETURN(0);
}
/*----------------------------------------------------------------------
* fill node_list_t struct from 1D surface file
*----------------------------------------------------------------------
*/
int sdata_from_1D ( opts_t * opts, param_t * p, node_list_t * N )
{
MRI_IMAGE * im;
float * fim;
int c;
ENTRY("sdata_from_1D");
if ( !opts || !N || !p )
RETURN(-1);
if ( (im = mri_read_1D(opts->sdata_file_1D)) == NULL )
{
fprintf(stderr,"** failed to read 1D file '%s'\n", opts->sdata_file_1D);
RETURN(-2);
}
/* transpose list to node-major order, and lose temp image */
p->sdata_im = mri_transpose(im);
mri_free(im);
if ( p->sdata_im->nx < 2 || p->sdata_im->ny < 1 )
{
fprintf(stderr,"** bad (%d x %d) surf data 1D file '%s'\n",
p->sdata_im->ny, p->sdata_im->nx, opts->sdata_file_1D);
RETURN(-3);
}
N->ilen = p->sdata_im->ny;
if ( opts->debug > 0 )
fprintf(stderr,"++ read 1D surface file '%s' (nx = %d, ny = %d)\n",
opts->sdata_file_1D, p->sdata_im->nx, p->sdata_im->ny );
/* only allocate space ilist */
if ( (N->ilist = (int *)malloc(N->ilen*sizeof(int))) == NULL )
{
fprintf(stderr,"** sf1D a: failed to allocate %d ints\n", N->ilen);
RETURN(-1);
}
/* first set the node index values */
fim = MRI_FLOAT_PTR( p->sdata_im );
for ( c = 0; c < N->ilen; c++, fim += p->sdata_im->nx )
N->ilist[c] = (int)*fim; /* set node index */
RETURN(0);
}
/*----------------------------------------------------------------------
* read surfaces (much stolen from SUMA_suma.c - thanks Ziad!)
*
* return 0 on success
*----------------------------------------------------------------------
*/
int fill_SUMA_structs ( opts_t * opts, SUMA_SurfSpecFile * spec )
{
int debug = 0, rv;
ENTRY("fill_SUMA_structs");
if ( opts->debug > 2 )
debug = 1;
if ( debug )
fputs( "-- SUMA_Create_CommonFields()...\n", stderr );
if ( opts->spec_file == NULL )
RETURN(-1);
/* initialize common fields struct */
SUMAg_CF = SUMA_Create_CommonFields();
if ( SUMAg_CF == NULL )
{
fprintf( stderr, "** failed SUMA_Create_CommonFields(), exiting...\n" );
RETURN(-1);
}
/* for SUMA type notifications */
if ( opts->debug > 3 )
{
SUMAg_CF->MemTrace = 1;
if ( opts->debug > 4 )
SUMAg_CF->InOut_Notify = 1;
}
if ( debug )
fputs( "-- SUMA_Alloc_DisplayObject_Struct()...\n", stderr );
SUMAg_DOv = SUMA_Alloc_DisplayObject_Struct(SUMA_MAX_DISPLAYABLE_OBJECTS);
if ( debug )
fputs( "-- SUMA_Read_SpecFile()...\n", stderr );
if ( SUMA_Read_SpecFile( opts->spec_file, spec) == 0 )
{
fprintf( stderr, "** failed SUMA_Read_SpecFile(), exiting...\n" );
RETURN(-1);
}
if ( debug )
SUMA_ShowSpecStruct(spec, stderr, 1);
rv = SUMA_spec_select_surfs(spec, opts->snames, S2V_MAX_SURFS, opts->debug);
if ( rv < 1 )
{
if ( rv == 0 )
fprintf(stderr,"** no named surfaces found in spec file\n");
RETURN(-1);
}
if ( opts->debug > 0 )
SUMA_ShowSpecStruct(spec, stderr, 1);
if ( SUMA_spec_set_map_refs(spec, opts->debug) != 0 )
RETURN(-1);
/* make sure only group was read from spec file */
if ( spec->N_Groups != 1 )
{
fprintf( stderr,"** error: N_Groups <%d> must be 1 in spec file <%s>\n",
spec->N_Groups, opts->spec_file );
RETURN(-1);
}
if ( debug )
fputs( "-- SUMA_LoadSpec_eng()...\n", stderr );
/* actually load the surface(s) from the spec file */
if (SUMA_LoadSpec_eng(spec,SUMAg_DOv,&SUMAg_N_DOv,opts->sv_file,debug,
SUMAg_CF->DsetList) == 0) /* DsetList - 26 Mar 2004 [ziad] */
{
fprintf( stderr, "** error: failed SUMA_LoadSpec_eng(), exiting...\n" );
RETURN(-1);
}
if ( opts->debug > 1 )
fprintf(stderr, "++ %d surfaces loaded.\n", spec->N_Surfs );
RETURN(0);
}
/*----------------------------------------------------------------------
* init_options - fill opts struct, display help
*----------------------------------------------------------------------
*/
int init_options ( opts_t * opts, int argc, char * argv [] )
{
int ac, ind;
ENTRY("init_options");
if ( argc < 2 )
{
usage( PROG_NAME, S2V_USE_LONG );
RETURN(-1);
}
/* clear out the options structure, pointers get explicit NULL */
memset( opts, 0, sizeof( opts_t) );
opts->gpar_file = NULL; opts->oset_file = NULL;
opts->spec_file = NULL; opts->sv_file = NULL;
opts->surf_xyz_1D_file = NULL; opts->sdata_file_1D = NULL;
opts->sdata_file_niml = NULL; opts->cmask_cmd = NULL;
opts->data_expr = NULL; opts->map_str = NULL;
opts->datum_str = NULL; opts->f_index_str = NULL;
opts->snames[0] = NULL; opts->snames[1] = NULL;
opts->dnode = -1; /* init to something invalid */
opts->dvox = -1; /* init to something invalid */
for ( ac = 1; ac < argc; ac++ )
{
/* alphabetical... */
if ( ! strncmp(argv[ac], "-cmask", 6) )
{
if ( (ac+1) >= argc )
{
fputs( "option usage: -cmask COMMAND\n\n", stderr );
usage( PROG_NAME, S2V_USE_SHORT );
RETURN(-1);
}
opts->cmask_cmd = argv[++ac];
}
else if ( ! strncmp(argv[ac], "-data_expr", 8) )
{
if ( (ac+1) >= argc )
{
fputs( "option usage: -data_expr EXPR\n\n", stderr );
usage( PROG_NAME, S2V_USE_SHORT );
RETURN(-1);
}
opts->data_expr = argv[++ac];
}
else if ( ! strncmp(argv[ac], "-datum", 6) )
{
if ( (ac+1) >= argc )
{
fputs( "option usage: -datum DTYPE\n\n", stderr );
usage( PROG_NAME, S2V_USE_SHORT );
RETURN(-1);
}
opts->datum_str = argv[++ac];
}
else if ( ! strncmp(argv[ac], "-debug", 6) )
{
if ( (ac+1) >= argc )
{
fputs( "option usage: -debug LEVEL\n\n", stderr );
usage( PROG_NAME, S2V_USE_SHORT );
RETURN(-1);
}
opts->debug = atoi(argv[++ac]);
if ( opts->debug < 0 || opts->debug > S2V_DEBUG_MAX_LEV )
{
fprintf( stderr, "bad debug level <%d>, should be in [0,%d]\n",
opts->debug, S2V_DEBUG_MAX_LEV );
usage( PROG_NAME, S2V_USE_SHORT );
RETURN(-1);
}
}
else if ( ! strncmp(argv[ac], "-dnode", 6) )
{
if ( (ac+1) >= argc )
{
fputs( "option usage: -dnode DEBUG_NODE\n\n", stderr );
usage( PROG_NAME, S2V_USE_SHORT );
RETURN(-1);
}
opts->dnode = atoi(argv[++ac]);
}
else if ( ! strncmp(argv[ac], "-dvoxel", 5) )
{
if ( (ac+1) >= argc )
{
fputs( "option usage: -dvoxel DEBUG_VOXEL\n\n", stderr );
usage( PROG_NAME, S2V_USE_SHORT );
RETURN(-1);
}
opts->dvox = atoi(argv[++ac]);
}
else if ( ! strncmp(argv[ac], "-f_index", 7) )
{
if ( (ac+1) >= argc )
{
fputs( "option usage: -f_index INDEX_TYPE\n\n", stderr );
usage( PROG_NAME, S2V_USE_SHORT );
RETURN(-1);
}
opts->f_index_str = argv[++ac];
}
else if ( ! strncmp(argv[ac], "-f_p1_fr", 9) )
{
if ( (ac+1) >= argc )
{
fputs( "option usage: -f_p1_fr FRACTION\n\n", stderr );
usage( PROG_NAME, S2V_USE_SHORT );
RETURN(-1);
}
opts->f_p1_fr = atof(argv[++ac]);
}
else if ( ! strncmp(argv[ac], "-f_pn_fr", 9) )
{
if ( (ac+1) >= argc )
{
fputs( "option usage: -f_pn_fr FRACTION\n\n", stderr );
usage( PROG_NAME, S2V_USE_SHORT );
RETURN(-1);
}
opts->f_pn_fr = atof(argv[++ac]);
}
else if ( ! strncmp(argv[ac], "-f_p1_mm", 9) )
{
if ( (ac+1) >= argc )
{
fputs( "option usage: -f_p1_mm DISTANCE\n\n", stderr );
usage( PROG_NAME, S2V_USE_SHORT );
RETURN(-1);
}
opts->f_p1_mm = atof(argv[++ac]);
}
else if ( ! strncmp(argv[ac], "-f_pn_mm", 9) )
{
if ( (ac+1) >= argc )
{
fputs( "option usage: -f_pn_mm DISTANCE\n\n", stderr );
usage( PROG_NAME, S2V_USE_SHORT );
RETURN(-1);
}
opts->f_pn_mm = atof(argv[++ac]);
}
else if ( ! strncmp(argv[ac], "-f_steps", 6) )
{
if ( (ac+1) >= argc )
{
fputs( "option usage: -f_steps NUM_STEPS\n\n", stderr );
usage( PROG_NAME, S2V_USE_SHORT );
RETURN(-1);
}
opts->f_steps = atoi(argv[++ac]);
}
else if ( ! strncmp(argv[ac], "-grid_parent", 5)||
! strncmp(argv[ac], "-inset", 6) ||
! strncmp(argv[ac], "-input", 6) )
{
if ( (ac+1) >= argc )
{
fputs( "option usage: -grid_parent INPUT_DSET\n\n", stderr );
usage( PROG_NAME, S2V_USE_SHORT );
RETURN(-1);
}
opts->gpar_file = argv[++ac];
}
else if ( ! strncmp(argv[ac], "-help", 5) )
{
usage( PROG_NAME, S2V_USE_LONG );
RETURN(-1);
}
else if ( ! strncmp(argv[ac], "-hist", 5) )
{
usage( PROG_NAME, S2V_USE_HIST );
RETURN(-1);
}
else if ( ! strncmp(argv[ac], "-map_func", 4) ) /* mapping function */
{
if ( (ac+1) >= argc )
{
fputs( "option usage: -map_func FUNCTION\n\n", stderr );
RETURN(-1);
}
opts->map_str = argv[++ac]; /* store user string for now */
}
else if ( ! strncmp(argv[ac], "-noscale", 4) )
{
opts->noscale = 1;
}
else if ( ! strncmp(argv[ac], "-prefix", 4) )
{
if ( (ac+1) >= argc )
{
fputs( "option usage: -prefix OUTPUT_PREFIX\n\n", stderr );
usage( PROG_NAME, S2V_USE_SHORT );
RETURN(-1);
}
opts->oset_file = argv[++ac];
}
else if ( ! strncmp(argv[ac], "-sdata_1D", 9) )
{
if ( (ac+1) >= argc )
{
fputs( "option usage: -sdata_1D SURF_DATA.1D\n\n", stderr );
usage( PROG_NAME, S2V_USE_SHORT );
RETURN(-1);
}
opts->sdata_file_1D = argv[++ac];
}
else if ( ! strncmp(argv[ac], "-sdata_niml", 9) )
{
if ( (ac+1) >= argc )
{
fputs( "option usage: -sdata_niml SURF_DATA.niml\n\n", stderr );
usage( PROG_NAME, S2V_USE_SHORT );
RETURN(-1);
}
opts->sdata_file_niml = argv[++ac];
}
else if ( ! strncmp(argv[ac], "-spec", 5) )
{
if ( (ac+1) >= argc )
{
fputs( "option usage: -spec SPEC_FILE\n\n", stderr );
usage( PROG_NAME, S2V_USE_SHORT );
RETURN(-1);
}
opts->spec_file = argv[++ac];
}
else if ( ! strncmp(argv[ac], "-surf_xyz_1D", 9) )
{
if ( (ac+1) >= argc )
{
fputs( "option usage: -surf_xyz_1D NODE_FILE\n\n", stderr );
usage( PROG_NAME, S2V_USE_SHORT );
RETURN(-1);
}
opts->surf_xyz_1D_file = argv[++ac];
}
else if ( ! strncmp(argv[ac], "-surf_", 6) )
{
if ( (ac+1) >= argc )
{
fputs( "option usage: -surf_X SURF_NAME\n\n", stderr );
usage( PROG_NAME, S2V_USE_SHORT );
RETURN(-1);
}
ind = argv[ac][6] - 'A';
if ( (ind < 0) || (ind >= S2V_MAX_SURFS) )
{
fprintf(stderr,"** -surf_X option: '%s' out of range,\n"
" use one of '-surf_A' through '-surf_%c'\n",
argv[ac], 'A'+S2V_MAX_SURFS-1);
RETURN(-1);
}
opts->snames[ind] = argv[++ac];
}
else if ( ! strncmp(argv[ac], "-sv", 3) )
{
if ( (ac+1) >= argc )
{
fputs( "option usage: -sv SURFACE_VOLUME\n\n", stderr );
usage( PROG_NAME, S2V_USE_SHORT );
RETURN(-1);
}
opts->sv_file = argv[++ac];
}
else if ( ! strncmp(argv[ac], "-sxyz_orient_as_gpar", 20) )
{
opts->sxyz_ori_gpar = 1;
}
else if ( ! strncmp(argv[ac], "-version", 2) )
{
usage( PROG_NAME, S2V_USE_VERSION );
RETURN(-1);
}
else /* invalid option */
{
fprintf( stderr, "invalid option <%s>\n", argv[ac] );
usage( PROG_NAME, S2V_USE_SHORT );
RETURN(-1);
}
}
if ( opts->debug > 1 )
disp_opts_t ( "++ opts read: ", opts );
RETURN(0);
}
/*----------------------------------------------------------------------
* validate_options - fill param struct from options
*
* - validate datasets
* - validate surface
*----------------------------------------------------------------------
*/
int validate_options ( opts_t * opts, param_t * p )
{
ENTRY("validate_options");
/* clear param struct - pointer get explicit NULL */
memset( p, 0, sizeof(*p) );
p->gpar = NULL; p->oset = NULL;
p->sxyz_im = NULL; p->sdata_im = NULL;
p->parser.pcode = NULL; p->cmask = NULL;
if ( check_map_func( opts->map_str ) == E_SMAP_INVALID )
RETURN(-1);
if ( validate_datasets( opts, p ) != 0 )
RETURN(-1);
if ( validate_surface( opts, p ) != 0 )
RETURN(-1);
if ( opts->data_expr )
{
PARSER_set_printout(1);
p->parser.pcode = PARSER_generate_code( opts->data_expr );
if ( p->parser.pcode == NULL )
{
fprintf(stderr,"** failed to generate code from expression '%s'\n",
opts->data_expr);
RETURN(-1);
}
PARSER_mark_symbols( p->parser.pcode, p->parser.has_sym );
if ( opts->debug > 1 )
disp_parser_t( "-- PARSER expr okay: ", &p->parser );
}
if ( opts->debug > 1 )
disp_param_t( "++ params set: ", p );
RETURN(0);
}
/*----------------------------------------------------------------------
* check_datum_type - determine type for output dataset
*
* currently allowable types are MRI_byte, MRI_short, MRI_float, MRI_double
*
* return datum type : on success
* -1 : on failure (no valid type can be determined)
*----------------------------------------------------------------------
*/
int check_datum_type ( char * datum_str, int default_type )
{
int c, dt = -1;
ENTRY("check_datum_type");
if ( datum_str )
{
for ( c = 0; c <= MRI_rgba; c++ )
if ( ! strcmp( datum_str, MRI_TYPE_name[c] ) )
{
dt = c;
break;
}
/* if we didn't find the requested type, inform the user and quit */
if ( c > MRI_rgba )
{
fprintf( stderr, "** invalid datum type name '%s'\n",
datum_str );
RETURN(-1);
}
}
else
{
dt = default_type;
}
if ( ( dt != MRI_byte ) &&
( dt != MRI_short ) &&
( dt != MRI_float ) &&
( dt != MRI_double ) )
{
fprintf( stderr, "** data type '%s' is not supported\n",
MRI_TYPE_name[dt] );
RETURN(-1);
}
RETURN(dt);
}
/*----------------------------------------------------------------------
* check_map_func
*
* - check for map_str
* - validate the map type
*----------------------------------------------------------------------
*/
int check_map_func ( char * map_str )
{
int map;
ENTRY("check_map_func");
if ( map_str == NULL )
{
fprintf( stderr, "** missing option: '-map_func FUNCTION'\n" );
RETURN(E_SMAP_INVALID);
}
map = s2v_map_type( map_str );
if ( map == E_SMAP_INVALID )
{
fprintf( stderr, "** invalid map string '%s'\n", map_str );
RETURN(-1);
}
RETURN(map);
}
/*----------------------------------------------------------------------
* s2v_map_type - return an E_SMAP_XXX code
*
* on failure, return -1 (E_SMAP_INVALID)
* else return >0 (a valid map code)
*----------------------------------------------------------------------
*/
int s2v_map_type ( char * map_str )
{
s2v_map_num map;
ENTRY("s2v_map_type");
if ( map_str == NULL )
{
fprintf( stderr, "** s2v_map_type: missing map_str parameter\n" );
RETURN((int)E_SMAP_INVALID);
}
if ( sizeof(gs2v_map_names) / sizeof(char *) != (int)E_SMAP_FINAL )
{
fprintf( stderr, "** error: gs2v_map_names/s2v_map_num mis-match\n");
RETURN((int)E_SMAP_INVALID);
}
for ( map = E_SMAP_MASK; map < E_SMAP_FINAL; map++ )
if ( !strcmp( map_str, gs2v_map_names[map] ) )
RETURN((int)map);
RETURN((int)E_SMAP_INVALID);
}
/*----------------------------------------------------------------------
* validate_surface
*
* Verify that the user entered options for both the spec file and
* the surface volume (AFNI dataset).
*----------------------------------------------------------------------
*/
int validate_surface ( opts_t * opts, param_t * p )
{
int errs = 0;
ENTRY("validate_surface");
if ( ! opts->surf_xyz_1D_file ) /* then check the surface input */
{
if ( opts->spec_file == NULL )
{
fprintf( stderr, "** missing '-spec_file SPEC_FILE' option\n" );
errs++;
}
if ( opts->sv_file == NULL )
{
fprintf( stderr, "** missing '-sv SURF_VOL' option\n" );
errs++;
}
if ( opts->snames[0] == NULL )
{
fprintf( stderr, "** missing '-surf_A SURF_NAME' option\n" );
errs++;
}
}
else if ( opts->spec_file != NULL )
{
fprintf(stderr,
"** cannot use both spec and xyz surface files, %s and %s\n",
opts->spec_file, opts->surf_xyz_1D_file);
errs++;
}
if ( opts->sdata_file_1D && opts->sdata_file_niml )
{
fprintf(stderr,"** cannot use both NIML and 1D surface files\n");
errs++;
}
/* rcr - hopefully this will disappear someday */
if ( opts->sdata_file_niml )
{
fprintf(stderr,"** sorry, the niml feature is coming soon...\n");
errs++;
}
if ( errs > 0 )
RETURN(-1);
RETURN(0);
}
/*----------------------------------------------------------------------
* validate_datasets
*
* Note that we do not validate the SURFACE_VOLUME AFNI dataset here.
* That is done in SUMA_LoadSpec_eng().
*
* Verify the AFNI dataset used for value output.
* Check for a cmask dataset and command.
* Verify that AFNI dataset and the mask have the same size.
*----------------------------------------------------------------------
*/
int validate_datasets( opts_t * opts, param_t * p )
{
ENTRY("validate_datasets");
p->gpar = THD_open_dataset( opts->gpar_file );
if ( !ISVALID_DSET(p->gpar) )
{
if ( opts->gpar_file == NULL )
fprintf( stderr, "** error: missing '-grid_parent DSET' option\n" );
else
fprintf( stderr, "** error: invalid input dataset '%s'\n",
opts->gpar_file);
RETURN(-1);
}
else if ( DSET_BRICK_TYPE(p->gpar, 0) == MRI_complex )
{
fprintf(stderr,
"** failure: cannot deal with complex-valued dataset, '%s'\n",
opts->gpar_file);
RETURN(-1);
}
p->nvox = DSET_NVOX( p->gpar );
set_3dmm_bounds( p->gpar, &p->f3mm_min, &p->f3mm_max );
if ( ! THD_filename_ok( opts->oset_file ) )
{
fprintf( stderr, "** illegal output prefix: '%s'\n",
opts->oset_file ? opts->oset_file : "<none>" );
RETURN(-1);
}
/* -------------------------------------------------------------------- */
/* check for cmask - casually stolen from 3dmaskdump.c (thanks, Bob! :) */
if ( opts->cmask_cmd != NULL )
{
int clen = strlen( opts->cmask_cmd );
char * cmd;
/* save original cmask command, as EDT_calcmask() is destructive */
cmd = (char *)malloc((clen + 1) * sizeof(char));
strcpy( cmd, opts->cmask_cmd );
p->cmask = EDT_calcmask( cmd, &p->ncmask );
free( cmd ); /* free EDT_calcmask() string */
if ( p->cmask == NULL )
{
fprintf( stderr, "** failure: cannot compute mask from option:\n"
" -cmask '%s'\n", opts->cmask_cmd );
RETURN(-1);
}
if ( p->ncmask != p->nvox )
{
fprintf( stderr, "** error: input and cmask datasets do not have "
"the same dimensions\n" );
RETURN(-1);
}
if ( ( p->ccount = THD_countmask( p->ncmask, p->cmask ) ) <= 0 )
{
fprintf( stderr, "** Warning! No voxels in computed cmask!\n" );
/* return -1; continue, and let the user deal with it... */
}
}
if ( opts->debug > 0 )
{
fprintf( stderr, "++ input dset has nvox = %d, nvals = %d",
p->nvox, DSET_NVALS(p->gpar) );
if ( p->cmask == NULL )
fputc( '\n', stderr );
else
fprintf( stderr, " (%d voxels in mask)\n", p->ccount );
}
RETURN(0);
}
/*----------------------------------------------------------------------
* usage - output usage information
*
* S2V_USE_SHORT - display brief output
* S2V_USE_LONG - display long output
* S2V_USE_VERSION - show the VERSION of the program
*----------------------------------------------------------------------
*/
int usage ( char * prog, int level )
{
ENTRY("usage");
if ( level == S2V_USE_SHORT )
fprintf(stderr,
" usage: %s [options] -spec SPEC_FILE -surf_A SURF_NAME \\\n"
" -grid_parent AFNI_DSET -sv SURF_VOL \\\n"
" -map_func MAP_FUNC -prefix OUTPUT_DSET\n"
"usage: %s -help\n",
prog, prog );
else if ( level == S2V_USE_LONG )
{
printf(
"\n"
"%s - map data from a surface domain to an AFNI volume domain\n"
"\n"
" usage: %s [options] -spec SPEC_FILE -surf_A SURF_NAME \\\n"
" -grid_parent AFNI_DSET -sv SURF_VOL \\\n"
" -map_func MAP_FUNC -prefix OUTPUT_DSET\n"
"\n"
" This program is meant to take as input a pair of surfaces,\n"
" optionally including surface data, and an AFNI grid parent\n"
" dataset, and to output a new AFNI dataset consisting of the\n"
" surface data mapped to the dataset grid space. The mapping\n"
" function determines how to map the surface values from many\n"
" nodes to a single voxel.\n"
"\n"
" Surfaces (from the spec file) are specified using '-surf_A'\n"
" (and '-surf_B', if a second surface is input). If two\n"
" surfaces are input, then the computed segments over node\n"
" pairs will be in the direction from surface A to surface B.\n"
"\n"
" The basic form of the algorithm is:\n"
"\n"
" o for each node pair (or single node)\n"
" o form a segment based on the xyz node coordinates,\n"
" adjusted by any '-f_pX_XX' options\n"
" o divide the segment up into N steps, according to \n"
" the '-f_steps' option\n"
" o for each segment point\n"
" o if the point is outside the space of the output\n"
" dataset, skip it\n"
" o locate the voxel in the output dataset which\n"
" corresponds to this segment point\n"
" o if the '-cmask' option was given, and the voxel\n"
" is outside the implied mask, skip it\n"
" o if the '-f_index' option is by voxel, and this\n"
" voxel has already been considered, skip it\n"
" o insert the surface node value, according to the\n"
" user-specified '-map_func' option\n"
"\n"
" Surface Coordinates:\n"
"\n"
" Surface coordinates are assumed to be in the Dicom\n"
" orientation. This information may come from the option\n"
" pair of '-spec' and '-sv', with which the user provides\n"
" the name of the SPEC FILE and the SURFACE VOLUME, along\n"
" with '-surf_A' and optionally '-surf_B', used to specify\n"
" actual surfaces by name. Alternatively, the surface\n"
" coordinates may come from the '-surf_xyz_1D' option.\n"
" See these option descriptions below.\n"
"\n"
" Note that the user must provide either the three options\n"
" '-spec', '-sv' and '-surf_A', or the single option,\n"
" '-surf_xyz_1D'.\n"
"\n"
" Surface Data:\n"
"\n"
" Surface domain data can be input via the '-sdata_1D'\n"
" option. In such a case, the data is with respect to the\n"
" input surface. The first column of the sdata_1D file\n"
" should be a node index, and following columns are that\n"
" node's data. See the '-sdata_1D' option for more info.\n"
"\n"
" If the surfaces have V values per node (pair), then the\n"
" resulting AFNI dataset will have V sub-bricks (unless the\n"
" user applies the '-data_expr' option).\n"
"\n"
" Mapping Functions:\n"
"\n"
" Mapping functions exist because a single volume voxel may\n"
" be occupied by multiple surface nodes or segment points.\n"
" Depending on how dense the surface mesh is, the number of\n"
" steps provided by the '-f_steps' option, and the indexing\n"
" type from '-f_index', even a voxel which is only 1 cubic\n"
" mm in volume may have quite a few contributing points.\n"
"\n"
" The mapping function defines how multiple surface values\n"
" are combined to get a single result in each voxel. For\n"
" example, the 'max' function will take the maximum of all\n"
" surface values contributing to each given voxel.\n"
"\n"
" Current mapping functions are listed under the '-map_func'\n"
" option, below.\n"
"\n"
"------------------------------------------------------------\n"
"\n"
" examples:\n"
"\n"
" 1. Map a single surface to an anatomical volume domain,\n"
" creating a simple mask of the surface. The output\n"
" dataset will be fred_surf+orig, and the orientation and\n"
" grid spacing will follow that of the grid parent. The\n"
" output voxels will be 1 where the surface exists, and 0\n"
" elsewhere.\n"
"\n"
" %s \\\n"
" -spec fred.spec \\\n"
" -surf_A pial \\\n"
" -sv fred_anat+orig \\\n"
" -grid_parent fred_anat+orig \\\n"
" -map_func mask \\\n"
" -prefix fred_surf\n"
"\n"
" 2. Map the cortical grey ribbon (between the white matter\n"
" surface and the pial surface) to an AFNI volume, where\n"
" the resulting volume is restriced to the mask implied by\n"
" the -cmask option.\n"
"\n"
" Surface data will come from the file sdata_10.1D, which\n"
" has 10 values per node, and lists only a portion of the\n"
" entire set of surface nodes. Each node pair will be form\n"
" a segment of 15 equally spaced points, the values from\n"
" which will be applied to the output dataset according to\n"
" the 'ave' filter. Since the index is over points, each\n"
" of the 15 points will have its value applied to the\n"
" appropriate voxel, even multiple times. This weights the\n"
" resulting average by the fraction of each segment that\n"
" occupies a given voxel.\n"
"\n"
" The output dataset will have 10 sub-bricks, according to\n"
" the 10 values per node index in sdata_10.1D.\n"
"\n"
" %s \\\n"
" -spec fred.spec \\\n"
" -surf_A smoothwm \\\n"
" -surf_B pial \\\n"
" -sv fred_anat+orig \\\n"
" -grid_parent 'fred_func+orig[0]' \\\n"
" -cmask '-a fred_func+orig[2] -expr step(a-0.6)' \\\n"
" -sdata_1D sdata_10.1D \\\n"
" -map_func ave \\\n"
" -f_steps 15 \\\n"
" -f_index points \\\n"
" -prefix fred_surf_ave\n"
"\n"
" 3. The inputs in this example are identical to those in\n"
" example 2, including the surface dataset, sdata_10.1D.\n"
" Again, the output dataset will have 10 sub-bricks.\n"
"\n"
" The surface values will be applied via the 'max_abs'\n"
" filter, with the intention of assigning to each voxel the\n"
" node value with the most significance. Here, the index\n"
" method does not matter, so it is left as the default,\n"
" 'voxel'.\n"
"\n"
" In this example, each node pair segment will be extended\n"
" by 20%% into the white matter, and by 10%% outside of the\n"
" grey matter, generating a \"thicker\" result.\n"
"\n"
" %s \\\n"
" -spec fred.spec \\\n"
" -surf_A smoothwm \\\n"
" -surf_B pial \\\n"
" -sv fred_anat+orig \\\n"
" -grid_parent 'fred_func+orig[0]' \\\n"
" -cmask '-a fred_func+orig[2] -expr step(a-0.6)' \\\n"
" -sdata_1D sdata_10.1D \\\n"
" -map_func max_abs \\\n"
" -f_steps 15 \\\n"
" -f_p1_fr -0.2 \\\n"
" -f_pn_fr 0.1 \\\n"
" -prefix fred_surf_max_abs\n"
"\n"
" 4. This is simliar to example 2. Here, the surface nodes\n"
" (coordinates) come from 'surf_coords_2.1D'. But these\n"
" coordinates do not happen to be in Dicomm orientation,\n"
" they are in the same orientation as the grid parent, so\n"
" the '-sxyz_orient_as_gpar' option is applied.\n"
"\n"
" Even though the data comes from 'sdata_10.1D', the output\n"
" AFNI dataset will only have 1 sub-brick. That is because\n"
" of the '-data_expr' option. Here, each applied surface\n"
" value will be the average of the sines of the first 3\n"
" data values (columns of sdata_10.1D).\n"
"\n"
" %s \\\n"
" -surf_xyz_1D surf_coords_2.1D \\\n"
" -sxyz_orient_as_gpar \\\n"
" -grid_parent 'fred_func+orig[0]' \\\n"
" -sdata_1D sdata_10.1D \\\n"
" -data_expr '(sin(a)+sin(b)+sin(c))/3' \\\n"
" -map_func ave \\\n"
" -f_steps 15 \\\n"
" -f_index points \\\n"
" -prefix fred_surf_ave_sine\n"
"\n"
" 5. In this example, voxels will get the maximum value from\n"
" column 3 of sdata_10.1D (as usual, column 0 is used for\n"
" node indices). The output dataset will have 1 sub-brick.\n"
"\n"
" Here, the output dataset is forced to be of type 'short',\n"
" regardless of what the grid parent is. Also, there will\n"
" be no scaling factor applied.\n"
"\n"
" To track the numbers for surface node #1234, the '-dnode'\n"
" option has been used, along with '-debug'. Additionally,\n"
" '-dvoxel' is used to track the results for voxel #6789.\n"
"\n"
" %s \\\n"
" -spec fred.spec \\\n"
" -surf_A smoothwm \\\n"
" -surf_B pial \\\n"
" -sv fred_anat+orig \\\n"
" -grid_parent 'fred_func+orig[0]' \\\n"
" -sdata_1D sdata_10.1D'[0,3]' \\\n"
" -map_func max \\\n"
" -f_steps 15 \\\n"
" -datum short \\\n"
" -noscale \\\n"
" -debug 2 \\\n"
" -dnode 1234 \\\n"
" -dvoxel 6789 \\\n"
" -prefix fred_surf_max\n"
"\n"
"------------------------------------------------------------\n"
"\n"
" REQUIRED COMMAND ARGUMENTS:\n"
"\n"
" -spec SPEC_FILE : SUMA spec file\n"
"\n"
" e.g. -spec fred.spec\n"
"\n"
" The surface specification file contains the list of\n"
" mappable surfaces that are used.\n"
"\n"
" See @SUMA_Make_Spec_FS and @SUMA_Make_Spec_SF.\n"
"\n"
" Note: this option, along with '-sv', may be replaced\n"
" by the '-surf_xyz_1D' option.\n"
"\n"
" -surf_A SURF_NAME : specify surface A (from spec file)\n"
" -surf_B SURF_NAME : specify surface B (from spec file)\n"
"\n"
" e.g. -surf_A smoothwm\n"
" e.g. -surf_A lh.smoothwm\n"
" e.g. -surf_B lh.pial\n"
"\n"
" This parameter is used to tell the program with surfaces\n"
" to use. The '-surf_A' parameter is required, but the\n"
" '-surf_B' parameter is an option.\n"
"\n"
" The surface names must uniquely match those in the spec\n"
" file, though a sub-string match is good enough. The\n"
" surface names are compared with the names of the surface\n"
" node coordinate files.\n"
"\n"
" For instance, given a spec file that has only the left\n"
" hemisphere in it, 'pial' should produce a unique match\n"
" with lh.pial.asc. But if both hemispheres are included,\n"
" then 'pial' would not be unique (matching rh.pail.asc,\n"
" also). In that case, 'lh.pial' would be better.\n"
"\n"
" -sv SURFACE_VOLUME : AFNI dataset\n"
"\n"
" e.g. -sv fred_anat+orig\n"
"\n"
" This is the AFNI dataset that the surface is mapped to.\n"
" This dataset is used for the initial surface node to xyz\n"
" coordinate mapping, in the Dicom orientation.\n"
"\n"
" Note: this option, along with '-spec', may be replaced\n"
" by the '-surf_xyz_1D' option.\n"
"\n"
" -surf_xyz_1D SXYZ_NODE_FILE : 1D coordinate file\n"
"\n"
" e.g. -surf_xyz_1D my_surf_coords.1D\n"
"\n"
" This ascii file contains a list of xyz coordinates to be\n"
" considered as a surface, or 2 sets of xyz coordinates to\n"
" considered as a surface pair. As usual, these points\n"
" are assumed to be in Dicom orientation. Another option\n"
" for coordinate orientation is to use that of the grid\n"
" parent dataset. See '-sxyz_orient_as_gpar' for details.\n"
"\n"
" This option is an alternative to the pair of options, \n"
" '-spec' and '-sv'.\n"
"\n"
" The number of rows of the file should equal the number\n"
" of nodes on each surface. The number of columns should\n"
" be either 3 for a single surface, or 6 for two surfaces.\n"
" \n"
" sample line of an input file (one surface):\n"
" \n"
" 11.970287 2.850751 90.896111\n"
" \n"
" sample line of an input file (two surfaces):\n"
" \n"
" 11.97 2.85 90.90 12.97 2.63 91.45\n"
" \n"
"\n"
" -grid_parent AFNI_DSET : AFNI dataset\n"
"\n"
" e.g. -grid_parent fred_function+orig\n"
"\n"
" This dataset is used as a grid and orientation master\n"
" for the output AFNI dataset.\n"
"\n"
" -map_func MAP_FUNC : surface to dataset function\n"
"\n"
" e.g. -map_func max\n"
" e.g. -map_func mask -f_steps 20\n"
"\n"
" This function applies to the case where multiple data\n"
" points get mapped to a single voxel, which is expected\n"
" since surfaces tend to have a much higher resolution\n"
" than AFNI volumes. In the general case data points come\n"
" from each point on each partitioned line segment, with\n"
" one segment per node pair. Note that these segments may\n"
" have length zero, such as when only a single surface is\n"
" input.\n"
"\n"
" See \"Mapping Functions\" above, for more information.\n"
"\n"
" The current mapping function for one surface is:\n"
"\n"
" mask : For each xyz location, set the corresponding\n"
" voxel to 1.\n"
"\n"
" The current mapping functions for two surfaces are as\n"
" follows. These descriptions are per output voxel, and\n"
" over the values of all points mapped to a given voxel.\n"
"\n"
" mask2 : if any points are mapped to the voxel, set\n"
" the voxel value to 1\n"
"\n"
" ave : average all values\n"
"\n"
" count : count the number of mapped data points\n"
"\n"
" min : find the minimum value from all mapped points\n"
"\n"
" max : find the maximum value from all mapped points\n"
"\n"
" max_abs: find the number with maximum absolute value\n"
" (the resulting value will retain its sign)\n"
"\n"
" -prefix OUTPUT_PREFIX : prefix for the output dataset\n"
"\n"
" e.g. -prefix anat_surf_mask\n"
"\n"
" This is used to specify the prefix of the resulting AFNI\n"
" dataset.\n"
"\n"
" ------------------------------\n"
" SUB-SURFACE DATA FILE OPTIONS:\n"
"\n"
" -sdata_1D SURF_DATA.1D : 1D sub-surface file, with data\n"
"\n"
" e.g. -sdata_1D roi3.1D\n"
"\n"
" This is used to specify a 1D file, which contains\n"
" surface indices and data. The indices refer to the\n"
" surface(s) read from the spec file.\n"
" \n"
" The format of this data file is a surface index and a\n"
" list of data values on each row. To be a valid 1D file,\n"
" each row must have the same number of columns.\n"
"\n"
" ------------------------------\n"
" OPTIONS SPECIFIC TO SEGMENT SELECTION:\n"
"\n"
" (see \"The basic form of the algorithm\" for more details)\n"
"\n"
" -f_steps NUM_STEPS : partition segments\n"
"\n"
" e.g. -f_steps 10\n"
" default: -f_steps 2 (or 1, the number of surfaces)\n"
"\n"
" This option specifies the number of points to divide\n"
" each line segment into, before mapping the points to the\n"
" AFNI volume domain. The default is the number of input\n"
" surfaces (usually, 2). The default operation is to have\n"
" the segment endpoints be the actual surface nodes,\n"
" unless they are altered with the -f_pX_XX options.\n"
"\n"
" -f_index TYPE : index by points or voxels\n"
"\n"
" e.g. -f_index points\n"
" e.g. -f_index voxels\n"
" default: -f_index voxels\n"
"\n"
" Along a single segment, the default operation is to\n"
" apply only those points mapping to a new voxel. The\n"
" effect of the default is that a given voxel will have\n"
" at most one value applied per voxel pair.\n"
"\n"
" If the user applies this option with 'points' or 'nodes'\n"
" as the argument, then every point along the segment will\n"
" be applied. This may be preferred if, for example, the\n"
" user wishes to have the average weighted by the number\n"
" of points occupying a voxel, not just the number of node\n"
" pair segments.\n"
"\n"
" Note: the following -f_pX_XX options are used to alter the\n"
" locations of the segment endpoints, per node pair.\n"
" The segments are directed, from the node on the first\n"
" surface to the node on the second surface. To modify\n"
" the first endpoint, use a -f_p1_XX option, and use\n"
" -f_pn_XX to modify the second.\n"
"\n"
" -f_p1_fr FRACTION : offset p1 by a length fraction\n"
"\n"
" e.g. -f_p1_fr -0.2\n"
" e.g. -f_p1_fr -0.2 -f_pn_fr 0.2\n"
"\n"
" This option moves the first endpoint, p1, by a distance\n"
" of the FRACTION times the original segment length. If\n"
" the FRACTION is positive, it moves in the direction of\n"
" the second endpoint, pn.\n"
"\n"
" In the example, p1 is moved by 20%% away from pn, which\n"
" will increase the length of each segment.\n"
"\n"
" -f_pn_fr FRACTION : offset pn by a length fraction\n"
"\n"
" e.g. -f_pn_fr 0.2\n"
" e.g. -f_p1_fr -0.2 -f_pn_fr 0.2\n"
"\n"
" This option moves pn by a distance of the FRACTION times\n"
" the original segment length, in the direction from p1 to\n"
" pn. So a positive fraction extends the segment, and a\n"
" negative fraction reduces it.\n"
"\n"
" In the example above, using 0.2 adds 20%% to the segment\n"
" length past the original pn.\n"
"\n"
" -f_p1_mm DISTANCE : offset p1 by a distance in mm.\n"
"\n"
" e.g. -f_p1_mm -1.0\n"
" e.g. -f_p1_mm -1.0 -f_pn_fr 1.0\n"
"\n"
" This option moves p1 by DISTANCE mm., in the direction\n"
" of pn. If the DISTANCE is positive, the segment gets\n"
" shorter. If DISTANCE is negative, the segment will get\n"
" longer.\n"
"\n"
" In the example, p1 is moved away from pn, extending the\n"
" segment by 1 millimeter.\n"
"\n"
" -f_pn_mm DISTANCE : offset pn by a distance in mm.\n"
"\n"
" e.g. -f_pn_mm 1.0\n"
" e.g. -f_p1_mm -1.0 -f_pn_fr 1.0\n"
"\n"
" This option moves pn by DISTANCE mm., in the direction\n"
" from the first point to the second. So if DISTANCE is\n"
" positive, the segment will get longer. If DISTANCE is\n"
" negative, the segment will get shorter.\n"
"\n"
" In the example, pn is moved 1 millimeter farther from\n"
" p1, extending the segment by that distance.\n"
"\n"
" ------------------------------\n"
" GENERAL OPTIONS:\n"
"\n"
" -cmask MASK_COMMAND : command for dataset mask\n"
"\n"
" e.g. -cmask '-a fred_func+orig[2] -expr step(a-0.8)'\n"
"\n"
" This option will produce a mask to be applied to the\n"
" output dataset. Note that this mask should form a\n"
" single sub-brick.\n"
"\n"
" This option follows the style of 3dmaskdump (since the\n"
" code for it was, uh, borrowed from there (thanks Bob!)).\n"
"\n"
" See '3dmaskdump -help' for more information.\n"
"\n"
" -data_expr EXPRESSION : apply expression to surface input\n"
"\n"
" e.g. -data_expr 17\n"
" e.g. -data_expr '(a+b+c+d)/4'\n"
" e.g. -data_expr '(sin(a)+sin(b))/2'\n"
"\n"
" This expression is applied to the list of data values\n"
" from the surface data file input via '-sdata_1D'. The\n"
" expression is applied for each node or node pair, to the\n"
" list of data values corresponding to that node.\n"
"\n"
" The letters 'a' through 'z' may be used as input, and\n"
" refer to columns 1 through 26 of the data file (where\n"
" column 0 is a surface node index). The data file must\n"
" have enough columns to support the expression. Is is\n"
" valid to have a constant expression without a data file.\n"
"\n"
" -datum DTYPE : set data type in output dataset\n"
"\n"
" e.g. -datum short\n"
" default: same as that of grid parent\n"
"\n"
" This option specifies the data type for the output AFNI\n"
" dataset. Valid choices are byte, short and float, which\n"
" are 1, 2 and 4 bytes for each data point, respectively.\n"
"\n"
" -debug LEVEL : verbose output\n"
"\n"
" e.g. -debug 2\n"
"\n"
" This option is used to print out status information \n"
" during the execution of the program. Current levels are\n"
" from 0 to 5.\n"
"\n"
" -dnode DEBUG_NODE : extra output for that node\n"
"\n"
" e.g. -dnode 123456\n"
"\n"
" This option requests additional debug output for the\n"
" given surface node. This index is with respect to the\n"
" input surface (included in the spec file, or through the\n"
" '-surf_xyz_1D' option).\n"
"\n"
" This will have no effect without the '-debug' option.\n"
"\n"
" -dvoxel DEBUG_VOXEL : extra output for that voxel\n"
"\n"
" e.g. -dvoxel 234567\n"
"\n"
" This option requests additional debug output for the\n"
" given volume voxel. This 1-D index is with respect to\n"
" the output AFNI dataset. One good way to find a voxel\n"
" index to supply is from output via the '-dnode' option.\n"
"\n"
" This will have no effect without the '-debug' option.\n"
"\n"
" -hist : show revision history\n"
"\n"
" Display module history over time.\n"
"\n"
" -help : show this help\n"
"\n"
" If you can't get help here, please get help somewhere.\n"
"\n"
" -noscale : no scale factor in output dataset\n"
"\n"
" If the output dataset is an integer type (byte, shorts\n"
" or ints), then the output dataset may end up with a\n"
" scale factor attached (see 3dcalc -help). With this\n"
" option, the output dataset will not be scaled.\n"
"\n"
" -sxyz_orient_as_gpar : assume gpar orientation for sxyz\n"
"\n"
" This option specifies that the surface coordinate points\n"
" in the '-surf_xyz_1D' option file have the orientation\n"
" of the grid parent dataset.\n"
"\n"
" When the '-surf_xyz_1D' option is applied the surface\n"
" coordinates are assumed to be in Dicom orientation, by\n"
" default. This '-sxyz_orient_as_gpar' option overrides\n"
" the Dicom default, specifying that the node coordinates\n"
" are in the same orientation as the grid parent dataset.\n"
"\n"
" See the '-surf_xyz_1D' option for more information.\n"
"\n"
" -version : show version information\n"
"\n"
" Show version and compile date.\n"
"\n"
"------------------------------------------------------------\n"
"\n"
" Author: R. Reynolds - %s\n"
"\n"
" (many thanks to Z. Saad and R.W. Cox)\n"
"\n",
prog, prog,
prog, prog, prog, prog, prog,
VERSION );
}
else if ( level == S2V_USE_HIST )
fputs(g_history, stdout);
else if ( level == S2V_USE_VERSION )
printf( "%s : %s, compile date: %s\n", prog, VERSION, __DATE__ );
else
fprintf( stderr, "usage called with illegal level <%d>\n", level );
RETURN(-1);
}
/*----------------------------------------------------------------------
* disp_opts_t - display the contents of the opts_t struct
*----------------------------------------------------------------------
*/
int disp_opts_t ( char * info, opts_t * opts )
{
ENTRY("disp_opts_t");
if ( info )
fputs( info, stderr );
if ( opts == NULL )
{
fprintf(stderr, "disp_opts_t: opts == NULL\n" );
RETURN(-1);
}
fprintf(stderr,
"options struct at %p :\n"
" gpar_file = %s\n"
" oset_file = %s\n"
" spec_file = %s\n"
" sv_file = %s\n"
" surf_xyz_1D_file = %s\n"
" sdata_file_1D = %s\n"
" sdata_file_niml = %s\n"
" cmask_cmd = %s\n"
" data_expr = %s\n"
" map_str = %s\n"
" datum_str = %s\n"
" f_index_str = %s\n"
" sxyz_ori_gpar = %d\n"
" debug, dnode, dvox = %d, %d, %d\n"
" noscale, f_steps = %d, %d\n"
" f_p1_fr, f_pn_fr = %f, %f\n"
" f_p1_mm, f_pn_mm = %f, %f\n"
, opts,
CHECK_NULL_STR(opts->gpar_file), CHECK_NULL_STR(opts->oset_file),
CHECK_NULL_STR(opts->spec_file), CHECK_NULL_STR(opts->sv_file),
CHECK_NULL_STR(opts->surf_xyz_1D_file),
CHECK_NULL_STR(opts->sdata_file_1D),
CHECK_NULL_STR(opts->sdata_file_niml),
CHECK_NULL_STR(opts->cmask_cmd), CHECK_NULL_STR(opts->data_expr),
CHECK_NULL_STR(opts->map_str), CHECK_NULL_STR(opts->datum_str),
CHECK_NULL_STR(opts->f_index_str), opts->sxyz_ori_gpar,
opts->debug, opts->dnode, opts->dvox, opts->noscale, opts->f_steps,
opts->f_p1_fr, opts->f_pn_fr, opts->f_p1_mm, opts->f_pn_mm
);
RETURN(0);
}
/*----------------------------------------------------------------------
* disp_param_t - display the contents of the param_t struct
*----------------------------------------------------------------------
*/
int disp_param_t ( char * info, param_t * p )
{
ENTRY("disp_param_t");
if ( info )
fputs( info, stderr );
if ( p == NULL )
{
fprintf(stderr,"disp_param_t: p == NULL\n");
RETURN(-1);
}
fprintf(stderr,
"param_t struct at %p :\n"
" gpar : vcheck = %p : %s\n"
" oset : vcheck = %p : %s\n"
" sxyz_im, sdata_im = %p, %p\n"
" f3mm_min (xyz) = (%f, %f, %f)\n"
" f3mm_max (xyz) = (%f, %f, %f)\n"
" nvox, nsubs = %d, %d\n"
" cmask = %p\n"
" ncmask, ccount = %d, %d\n"
, p,
p->gpar, ISVALID_DSET(p->gpar) ? "valid" : "invalid",
p->oset, ISVALID_DSET(p->oset) ? "valid" : "invalid",
p->sxyz_im, p->sdata_im,
p->f3mm_min.xyz[0], p->f3mm_min.xyz[1], p->f3mm_min.xyz[2],
p->f3mm_max.xyz[0], p->f3mm_max.xyz[1], p->f3mm_max.xyz[2],
p->nvox, p->nsubs, p->cmask, p->ncmask, p->ccount
);
RETURN(0);
}
/*----------------------------------------------------------------------
* disp_node_list_t - display the contents of the node_list_t struct
*----------------------------------------------------------------------
*/
int disp_node_list_t ( char * info, node_list_t * d )
{
ENTRY("disp_node_list_t");
if ( info )
fputs( info, stderr );
if ( d == NULL )
{
fprintf(stderr,"disp_node_list_t: d == NULL\n");
RETURN(-1);
}
fprintf(stderr,
"node_list_t struct at %p :\n"
" nodes = %p\n"
" depth, nnodes = %d, %d\n"
" fdata, ilist = %p, %p\n"
" ilen = %d\n"
, d,
d->nodes, d->depth, d->nnodes,
d->fdata, d->ilist, d->ilen
);
RETURN(0);
}
/*----------------------------------------------------------------------
* disp_s2v_opts_t - display the contents of the s2v_opts_t struct
*----------------------------------------------------------------------
*/
int disp_s2v_opts_t ( char * info, s2v_opts_t * sopt )
{
ENTRY("disp_s2v_opts_t");
if ( info )
fputs( info, stderr );
if ( sopt == NULL )
{
fprintf(stderr,"disp_s2v_opts_t: sopt == NULL\n" );
RETURN(-1);
}
fprintf(stderr,
"s2v_opts_t struct at %p :\n"
" map, datum, noscale = %d, %d, %d\n"
" debug, dnode, dvox = %d, %d, %d\n"
" sxyz_ori_gpar = %d\n"
" cmask = %p\n"
" f_steps, f_index = %d, %d\n"
" f_p1_fr, f_pn_fr = %f, %f\n"
" f_p1_mm, f_pn_mm = %f, %f\n"
, sopt,
sopt->map, sopt->datum, sopt->noscale,
sopt->debug, sopt->dnode, sopt->dvox, sopt->sxyz_ori_gpar,
sopt->cmask, sopt->f_steps, sopt->f_index,
sopt->f_p1_fr, sopt->f_pn_fr, sopt->f_p1_mm, sopt->f_pn_mm
);
RETURN(0);
}
/*----------------------------------------------------------------------
* disp_parser_t - display the contents of the parser_t struct
*----------------------------------------------------------------------
*/
int disp_parser_t ( char * info, parser_t * d )
{
int c;
ENTRY("disp_parser_t");
if ( info )
fputs( info, stderr );
if ( d == NULL )
{
fprintf(stderr,"disp_parser_t: d == NULL\n" );
RETURN(-1);
}
fprintf(stderr, "parser_t struct at %p :\n"
" pcode = %p\n"
" max_sym = %d\n",
d, d->pcode, d->max_sym );
if ( d->pcode )
{
fprintf(stderr, " num_code = %d\n"
" c_code = %s\n",
d->pcode->num_code, d->pcode->c_code );
fprintf(stderr, " has_sym =");
for ( c = 0; c < 26; c++ )
fprintf(stderr, " %d", d->has_sym[c]);
fputc('\n', stderr);
}
RETURN(0);
}
/*----------------------------------------------------------------------
* set_3dmm_bounds - note 3dmm bounding values - v1.2
*
* This is an outer bounding box, like FOV, not SLAB.
*----------------------------------------------------------------------
*/
int set_3dmm_bounds ( THD_3dim_dataset *dset, THD_fvec3 *min, THD_fvec3 *max)
{
float tmp;
int c;
ENTRY("set_3dmm_bounds");
if ( !dset || !min || !max )
{
fprintf(stderr, "** invalid params to set_3dmm_bounds: (%p,%p,%p)\n",
dset, min, max );
RETURN(-1);
}
/* get undirected bounds */
min->xyz[0] = DSET_XORG(dset) - 0.5 * DSET_DX(dset);
max->xyz[0] = min->xyz[0] + DSET_NX(dset) * DSET_DX(dset);
min->xyz[1] = DSET_YORG(dset) - 0.5 * DSET_DY(dset);
max->xyz[1] = min->xyz[1] + DSET_NY(dset) * DSET_DY(dset);
min->xyz[2] = DSET_ZORG(dset) - 0.5 * DSET_DZ(dset);
max->xyz[2] = min->xyz[2] + DSET_NZ(dset) * DSET_DZ(dset);
for ( c = 0; c < 3; c++ )
if ( min->xyz[c] > max->xyz[c] )
{
tmp = min->xyz[c];
min->xyz[c] = max->xyz[c];
max->xyz[c] = tmp;
}
RETURN(0);
}
/*----------------------------------------------------------------------
* dist_f3mm - return Euclidean distance between the points
* - v1.2
*----------------------------------------------------------------------
*/
float dist_f3mm( THD_fvec3 * p1, THD_fvec3 * p2 )
{
double d0, d1, d2;
ENTRY("dist_f3mm");
if ( p1 == NULL || p2 == NULL )
{
fprintf( stderr, "** dist_f3mm: invalid params (%p,%p)\n", p1, p2 );
RETURN(0.0);
}
d0 = p1->xyz[0] - p2->xyz[0];
d1 = p1->xyz[1] - p2->xyz[1];
d2 = p1->xyz[2] - p2->xyz[2];
RETURN(sqrt(d0*d0 + d1*d1 + d2*d2));
}
/*----------------------------------------------------------------------
* f3mm_out_of_bounds - check wether cp is between min and max
* - v1.2 [rickr]
*----------------------------------------------------------------------
*/
int f3mm_out_of_bounds( THD_fvec3 * cp, THD_fvec3 * min, THD_fvec3 * max )
{
int c;
ENTRY("f3mm_out_of_bounds");
if ( !cp || !min || !max )
RETURN(-1);
for ( c = 0; c < 3; c++ )
{
if ( ( cp->xyz[c] < min->xyz[c] ) ||
( cp->xyz[c] > max->xyz[c] ) )
RETURN(-1);
}
RETURN(0);
}
/*----------------------------------------------------------------------
* integral_doubles - are all double values integral?
*----------------------------------------------------------------------
*/
int integral_doubles( double * dp, int nvals )
{
ENTRY("integral_doubles");
while ( nvals > 0 )
{
if ( ((double)(int)*dp) != *dp )
RETURN(0);
dp++;
nvals--;
}
RETURN(1);
}
syntax highlighted by Code2HTML, v. 0.9.1