do { if (!(ptr))                                                  \
             { fprintf(stderr,"** failed to allocate %d '%s' elements\n", \
                        nel, name);                                       \
               exit(1);                                                   \
             } } while (0)

00034 {
00035     E_SM_INVALID,                       /* do not change INVALID or FINAL     */
00036     E_SM_ANG_NORMS,                     /* give angle diff between normals    */
00037     E_SM_ANG_NS_A,                      /* give norm A angle to segment       */
00038     E_SM_ANG_NS_B,                      /* give norm B angle to segment       */
00039     E_SM_COORD_A,                       /* give coordinates of first node     */
00040     E_SM_COORD_B,                       /* give coordinates of second node    */
00041     E_SM_N_AREA_A,                      /* give node area of first node       */
00042     E_SM_N_AREA_B,                      /* give node area of second node      */
00043     E_SM_N_AVEAREA_A,                   /* average area of included triangles */
00044     E_SM_N_AVEAREA_B,                   /* average area of included triangles */
00045     E_SM_NTRI,                          /* give number of included triangles  */
00046     E_SM_NODE_VOL,                      /* give volume for each node          */
00047     E_SM_NODES,                         /* print out node indices             */
00048     E_SM_NORM_A,                        /* give normal vector for first node  */
00049     E_SM_NORM_B,                        /* give normal vector for second node */
00050     E_SM_THICK,                         /* display between surface thickness  */
00051     E_SM_FINAL
00052 } smeasure_codes_e;

01298 {
01299 ENTRY("add_to_flist");
01300 
01301     if ( F->nused >= F->nalloc )
01302     {
01303         F->nalloc += ST_DEFAULT_FALLOC;
01304         F->names   = (char **)realloc(F->names, F->nalloc*sizeof(char *));
01305         F->codes   = (int   *)realloc(F->codes, F->nalloc*sizeof(int   ));
01306 
01307         if ( !F->names || !F->codes )
01308         {
01309             fprintf(stderr,"** failed to allocate for %d functs (%p,%p)\n",
01310                     F->nalloc, F->names, F->codes);
01311             RETURN(-1);
01312         }
01313     }
01314 
01315     F->names[F->nused] = fname;
01316     F->codes[F->nused] = check_func_name(fname);
01317 
01318     if ( F->codes[F->nused] == E_SM_INVALID )
01319     {
01320         fprintf(stderr,"** function '%s' is not valid\n",
01321                 CHECK_NULL_STR(fname));
01322         RETURN(-1);
01323     }
01324 
01325     F->nused++;                 /* we've got another good one */
01326 
01327     RETURN(0);
01328 }

01154 {
01155     SUMA_SurfaceObject * so;
01156     int                  count;
01157 
01158 ENTRY("all_mappable_surfs");
01159 
01160     p->S.slist = (SUMA_SurfaceObject **)malloc(SUMAg_N_DOv *
01161                                                sizeof(SUMA_SurfaceObject *));
01162     ALLOC_CHECK(p->S.slist,"SUMA_SurfaceObject *",SUMAg_N_DOv);
01163 
01164     p->S.salloc = SUMAg_N_DOv;          /* number allocated for              */
01165     p->S.nsurf  = 0;                    /* number of mappable surfaces found */
01166 
01167     for ( count = 0; count < SUMAg_N_DOv; count++ )
01168     {
01169         if ( ! SUMA_isSO(SUMAg_DOv[count]) )
01170             continue;
01171 
01172         so = (SUMA_SurfaceObject *)SUMAg_DOv[count].OP;
01173 
01174 /*      if ( ! SUMA_isINHmappable(so) )       -  deprecated  [v1.5] */
01175 
01176 #if 0   /* do not require this [v1.9] */
01177 
01178         if ( ! so->AnatCorrect )
01179         {
01180             if ( opts->debug )
01181                 fprintf(stderr,"** warning: surface '%s' is not mappable, "
01182                         "skipping...\n", so->Label ? so->Label : "<unnamed>");
01183             if ( opts->debug > 1 )
01184                 fprintf(stderr,"** consider adding the following to the "
01185                                "surface definition in the spec file:\n"
01186                                "       Anatomical = Y\n");
01187             continue;
01188         }
01189 #endif
01190 
01191         if ( opts->debug > 1 )
01192         {
01193             fprintf(stderr,"-------- surface #%d (%s) --------\n",
01194                     p->S.nsurf, so->Label ? so->Label : "<unnamed>");
01195             if ( opts->debug > 3 )
01196                 SUMA_Print_Surface_Object(so, stderr);
01197         }
01198 
01199         p->S.slist[p->S.nsurf] = so;
01200         p->S.nsurf++;
01201     }
01202 
01203     if ( p->S.nsurf > 0 )
01204     {
01205         p->S.nnodes = p->S.slist[0]->N_Node;
01206         p->S.nfaces = p->S.slist[0]->N_FaceSet;
01207     }
01208 
01209     if ( opts->debug )
01210         fprintf(stderr, "++ found %d mappable surfaces\n", p->S.nsurf);
01211 
01212     RETURN(0);
01213 }

02180 {
02181     int fnum;
02182 
02183 ENTRY("check_func_name");
02184 
02185     if ( !func )
02186         RETURN(E_SM_INVALID);
02187 
02188     /* just to be safe, let's verify the names and enum */
02189 
02190     if ( sizeof(g_sm_names)/sizeof(char *) != (int)E_SM_FINAL )
02191     {
02192         fprintf(stderr,"** error: g_sm_names mis-match\n");
02193         RETURN(E_SM_INVALID);
02194     }
02195 
02196     if ( sizeof(g_sm_desc)/sizeof(char *) != (int)E_SM_FINAL )
02197     {
02198         fprintf(stderr,"** error: g_sm_desc mis-match\n");
02199         RETURN(E_SM_INVALID);
02200     }
02201 
02202     for ( fnum = E_SM_INVALID; fnum < E_SM_FINAL; fnum++ )
02203         if ( !strcmp(func, g_sm_names[fnum]) )
02204             RETURN(fnum);
02205 
02206     RETURN(E_SM_INVALID);
02207 }

01054 {
01055     SUMA_SurfaceObject    * so0, *so1;
01056     double                  totalv;
01057     double                  v0, v1, v2;
01058     float                 * fvlist, * xyzn0, * xyzn1;
01059     float                   vmin, vmax;
01060     int                   * fp;
01061     int                     i0, i1, i2;
01062     int                     face, nnodes, test;
01063 
01064 ENTRY("compute_face_vols");
01065 
01066     so0 = p->S.slist[0];
01067     so1 = p->S.slist[1];
01068 
01069     if ( !so0 || !so1 )
01070     {
01071         fprintf(stderr,"** cfv: missing SurfaceObject (%p, %p)\n", so0, so1 );
01072         RETURN(-1);
01073     }
01074 
01075     fp    = so0->FaceSetList;  /* we need only 1, as they are the same */
01076     xyzn0 = so0->NodeList;
01077     xyzn1 = so1->NodeList;
01078     
01079     if ( !fp || !xyzn0 || !xyzn1 )
01080     {
01081         fprintf(stderr,"** cfv: missing face set list data (%p,%p,%p)\n",
01082                 fp, xyzn0, xyzn1);
01083         RETURN(-1);
01084     }
01085 
01086     fvlist = (float *)malloc(p->S.nfaces*sizeof(float));
01087     ALLOC_CHECK(fvlist, "float", p->S.nfaces);
01088     p->S.fvol = fvlist;
01089     
01090     vmin = 10000.0;  vmax = -1.0;
01091     nnodes = p->S.nnodes;
01092     totalv = 0.0;
01093     for ( face = 0; face < p->S.nfaces; face++ )
01094     {
01095         i0 = fp[0];  i1 = fp[1];  i2 = fp[2];
01096 
01097         if ( i0 < 0 || i1 < 0 || i2 < 0 ||
01098              i0 >= nnodes || i1 >= nnodes || i2 >= nnodes )
01099         {
01100             fprintf(stderr,"** cfv: face %d, index out of range [0,%d]\n"
01101                            "        indices are (%d,%d,%d)\n",
01102                            face, nnodes-1, i0, i1, i2);
01103             RETURN(-1);
01104         }
01105 
01106         test = (i0 == opts->dnode || i1 == opts->dnode || i2 == opts->dnode);
01107 
01108         i0 *= 3;  i1 *= 3;  i2 *= 3;    /* node index -> (float*) index */
01109 
01110         v0 = tetra_volume( xyzn0+i0, xyzn0+i1, xyzn0+i2, xyzn1+i0 );
01111         v1 = tetra_volume( xyzn0+i1, xyzn0+i2, xyzn1+i0, xyzn1+i1 );
01112         v2 = tetra_volume( xyzn0+i2, xyzn1+i0, xyzn1+i1, xyzn1+i2 );
01113 
01114         if ( v0 < 0.0 || v1 < 0.0 || v2 < 0.0 ||
01115              (opts->debug > 2 && test) )
01116         {
01117             fprintf(stderr,"** cfv: check volume: %f = %f + %f + %f, face %d\n",
01118                            v0+v1+v2, v0, v1, v2, face );
01119             fprintf(stderr,"        nodes %d, %d, %d\n", i0/3, i1/3, i2/3);
01120 
01121             if ( v0 < 0.0 || v1 < 0.0 || v2 < 0.0 )
01122                 RETURN(-1);
01123         }
01124 
01125         fvlist[face] = v0 + v1 + v2;
01126         totalv += fvlist[face];
01127         if ( fvlist[face] < vmin )  vmin = fvlist[face];
01128         if ( fvlist[face] > vmax )  vmax = fvlist[face];
01129 
01130         fp += 3;
01131     }
01132 
01133     if ( opts->debug > 0 )
01134     {
01135         fprintf(stderr,"++ total face volume = %f\n", totalv);
01136         if ( opts->debug > 1 )
01137         {
01138             fprintf(stderr,"++ volumes: faces 0, 1 (of %d), vols = %f, %f\n",
01139                     p->S.nfaces, fvlist[0], fvlist[1]);
01140             fprintf(stderr,"-- faces: vmin, vmax = %f, %f\n", vmin, vmax);
01141         }
01142     }
01143 
01144     RETURN(0);
01145 }

00932 {
00933     SUMA_SurfaceObject    * so;
00934     double                  sum;
00935     float                 * alist;
00936     int                   * flist;
00937     int                     node, c;
00938 
00939 ENTRY("compute_node_areas");
00940 
00941     if ( sindex < 0 || sindex >= p->S.nsurf || sindex >= 2 )
00942     {
00943         fprintf(stderr,"** cna: surf index <%d> is out of range\n",sindex);
00944         RETURN(-1);
00945     }
00946 
00947     so = p->S.slist[sindex];            /* just for ease of typing */
00948 
00949     if ( ! so->PolyArea )
00950     {
00951         fprintf(stderr,"** cna: no PolyArea to compute from\n");
00952         RETURN(-1);
00953     }
00954 
00955     alist = (float *)malloc(p->S.nnodes*sizeof(float));
00956     ALLOC_CHECK(alist, "float", p->S.nnodes);
00957     p->S.narea[sindex] = alist;
00958 
00959     for ( node = 0; node < p->S.nnodes; node++ )
00960     {
00961         flist = so->MF->NodeMemberOfFaceSet[node];
00962         sum = 0.0;
00963         for (c = 0; c < so->MF->N_Memb[node]; c++)
00964         {
00965             if ( flist[c] < 0 || flist[c] >= so->N_FaceSet )
00966             {
00967                 fprintf(stderr,"** cna: FaceSet mis-match flist,max = %d,%d\n",
00968                         flist[c], so->N_FaceSet);
00969                 free(p->S.narea[sindex]);
00970                 RETURN(-1);
00971             }
00972 
00973             sum += so->PolyArea[flist[c]];
00974         }
00975 
00976         alist[node] = sum/3.0;
00977 
00978         if ( node == opts->dnode )
00979             fprintf(stderr, "-- dnode %d: area = %s (%d faces, surf %s)\n",
00980                     node, MV_format_fval(alist[node]),
00981                     so->MF->N_Memb[node], CHECK_NULL_STR(so->Label));
00982     }
00983 
00984     RETURN(0);
00985 }

00996 {
00997     SUMA_SurfaceObject    * so;
00998     double                  sum;
00999     float                 * nvols;
01000     int                   * flist;
01001     int                     node, c;
01002 
01003 ENTRY("compute_node_vols");
01004 
01005     so = p->S.slist[0];                 /* just for ease of typing */
01006 
01007     nvols = (float *)malloc(p->S.nnodes*sizeof(float));
01008     ALLOC_CHECK(nvols, "float", p->S.nnodes);
01009     p->S.nvol = nvols;
01010 
01011     for ( node = 0; node < p->S.nnodes; node++ )
01012     {
01013         if ( node == opts->dnode && opts->debug > 1 )
01014             fprintf(stderr,"-- dnode %d, facelist (%d) :\n        ",
01015                     node, so->MF->N_Memb[node]);
01016         flist = so->MF->NodeMemberOfFaceSet[node];
01017         sum = 0.0;
01018         for (c = 0; c < so->MF->N_Memb[node]; c++)
01019         {
01020             if ( flist[c] < 0 || flist[c] >= so->N_FaceSet )
01021             {
01022                 fprintf(stderr,"** cnv: FaceSet mis-match flist,max = %d,%d\n",
01023                         flist[c], so->N_FaceSet);
01024                 free(p->S.nvol);
01025                 RETURN(-1);
01026             }
01027             if ( node == opts->dnode && opts->debug > 1 )
01028                 fprintf(stderr,"  %d", flist[c]);
01029 
01030             sum += p->S.fvol[flist[c]];
01031         }
01032 
01033         nvols[node] = sum/3.0;
01034 
01035         if ( node == opts->dnode && opts->debug > 0 )
01036             fprintf(stderr,"\n-- volume = %f\n", sum);
01037     }
01038 
01039     RETURN(0);
01040 }

00899 {
00900     res[0] = u[1]*v[2] - u[2]*v[1];
00901     res[1] = u[2]*v[0] - u[0]*v[2];
00902     res[2] = u[0]*v[1] - u[1]*v[0];
00903 
00904     return 0;
00905 }

02446 {
02447     if ( info )
02448         fputs( info, stderr );
02449 
02450     if ( ! d )
02451     {
02452         fprintf(stderr,"** disp_f3_point: d == NULL\n");
02453         return -1;
02454     }
02455 
02456     fprintf(stderr,"(%6s, ", MV_format_fval(d[0]));
02457     fprintf(stderr,"%6s, ",  MV_format_fval(d[1]));
02458     fprintf(stderr,"%6s)\n", MV_format_fval(d[2]));
02459 
02460     return 0;
02461 }

02352 {
02353     if ( info )
02354         fputs( info, stderr );
02355 
02356     if ( ! d )
02357     {
02358         fprintf(stderr,"** disp_func_t: d == NULL\n");
02359         return -1;
02360     }
02361 
02362     fprintf(stderr,
02363             "func_t struct at %p:\n"
02364             "    names, codes   = %p, %p\n"
02365             "    nalloc, nused  = %d, %d\n",
02366             d,
02367             d->names, d->codes, d->nalloc, d->nused);
02368 
02369     return 0;
02370 }

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 }

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 }

02319 {
02320     if ( info )
02321         fputs( info, stderr );
02322 
02323     if ( ! d )
02324     {
02325         fprintf(stderr,"** disp_surf_t: d == NULL\n");
02326         return -1;
02327     }
02328 
02329     fprintf(stderr,
02330             "surf_t struct at %p:\n"
02331             "    spec N_Surfs   = %d\n"
02332             "    spec N_Groups  = %d\n"
02333             "    slist          = %p\n"
02334             "    narea[0,1]     = %p, %p\n"
02335             "    nvol, fvol     = %p, %p\n"
02336             "    nsurf, salloc  = %d, %d\n"
02337             "    nnodes, nfaces = %d, %d\n",
02338             d,
02339             d->spec.N_Surfs, d->spec.N_Groups, d->slist,
02340             d->narea[0], d->narea[1], d->nvol, d->fvol,
02341             d->nsurf, d->salloc, d->nnodes, d->nfaces);
02342 
02343     return 0;
02344 }

02292 {
02293     double diff, sum;
02294     int    c;
02295 
02296     if ( len < 1 || p1 == NULL || p2 == NULL )
02297     {
02298         fprintf(stderr, "** dist_fn: invalid params (%d,%p,%p)\n",
02299                 len, p1, p2);
02300         return 0.0;
02301     }
02302 
02303     sum = 0.0;
02304     for ( c = 0; c < len; c++ )
02305     {
02306         diff = p1[c] - p2[c];
02307         sum += diff * diff;
02308     }
02309 
02310     return sqrt(sum);
02311 }

00913 {
00914     double res;
00915 
00916     res  = u[0]*v[0];
00917     res += u[1]*v[1];
00918     res += u[2]*v[2];
00919 
00920     return res;
00921 }

02215 {
02216 
02217 ENTRY("final_cleanup");
02218 
02219     /* first, close the output file, the rest are in order */
02220     if ( p->outfp != stdout )
02221         fclose(p->outfp);
02222 
02223     if ( p->S.narea[0] )  free(p->S.narea[0]);
02224     if ( p->S.narea[1] )  free(p->S.narea[1]);
02225     if ( p->S.slist    )  free(p->S.slist);
02226     if ( p->S.nvol     )  free(p->S.nvol);
02227     if ( p->S.fvol     )  free(p->S.fvol);
02228 
02229     if ( p->F->nalloc > 0 )
02230     {
02231         free(p->F->names);
02232         free(p->F->codes);
02233     }
02234 
02235     if ( p->nodes )
02236         free(p->nodes);
02237 
02238     if ( p->cmask )
02239         free(p->cmask);
02240 
02241     if ( opts->debug > 2 )
02242         fprintf(stderr,"-- freeing SUMA data...\n");
02243 
02244     if ( ( SUMAg_DOv != NULL ) &&
02245          ( SUMA_Free_Displayable_Object_Vect(SUMAg_DOv, SUMAg_N_DOv) == 0 ) )
02246         fprintf(stderr, "** failed SUMA_Free_Displayable_Object_Vect()\n" );
02247 
02248     if ( ( SUMAg_SVv != NULL ) &&
02249          ( SUMA_Free_SurfaceViewer_Struct_Vect(SUMAg_SVv, SUMAg_N_SVv) == 0 ) )
02250         fprintf( stderr, "** failed SUMA_Free_SurfaceViewer_Struct_Vect()\n" );
02251 
02252     if ( ( SUMAg_CF != NULL ) && ( SUMA_Free_CommonFields(SUMAg_CF) == 0 ) )
02253         fprintf( stderr, "** failed SUMA_Free_CommonFields()\n" );
02254 
02255     RETURN(0);
02256 }

00446 {
00447     double sums;
00448     int    c;
00449 
00450 ENTRY("fvec_magnitude");
00451 
00452     sums = 0.0;
00453     for ( c = 0; c < length; c++ )
00454         sums += v[c] * v[c];
00455 
00456     RETURN(sqrt(sums));
00457 }

01605 {
01606     char * cmd;
01607     int    clen;
01608 
01609 ENTRY("get_cmask");
01610 
01611     if ( ! opts->cmask_cmd )
01612         RETURN(0);
01613 
01614     clen = strlen(opts->cmask_cmd);
01615     cmd  = (char *)malloc((clen + 1)*sizeof(char));
01616     ALLOC_CHECK(cmd, "char", clen);
01617 
01618     strcpy(cmd, opts->cmask_cmd);
01619 
01620     p->cmask = EDT_calcmask(cmd, &p->ncmask);
01621 
01622     free(cmd);          /* we are done with the now corrupted command */
01623 
01624     if ( ! p->cmask || p->ncmask < 1 )
01625     {
01626         fprintf(stderr,"** failure: cannot compute cmask from option:\n"
01627                 "   -cmask '%s'\n", opts->cmask_cmd);
01628         RETURN(-1);
01629     }
01630 
01631     p->ccount = THD_countmask( p->ncmask, p->cmask );
01632 
01633     if ( p->ccount < 1 )                /* do not quit */
01634         fprintf(stderr,"** warning: cmask is empty from option\n"
01635                 "   -cmask '%s'\n", opts->cmask_cmd);
01636 
01637     if ( opts->debug > 0 )
01638         fprintf(stderr,"++ have cmask with %d of %d set entries\n",
01639                 p->ccount, p->ncmask);
01640 
01641     RETURN(0);
01642 }

00622 {
00623     int rv;
00624 
00625 ENTRY("get_surf_data");
00626 
00627     rv = spec2SUMA(&p->S.spec, opts);
00628     if ( rv != 0 )
00629         RETURN(rv);
00630 
00631     if ( opts->debug > 1 )
00632         fprintf(stderr,"-- timing: Surf (spec2SUMA)  : time = %.3f\n",
00633                 COX_clock_time());
00634 
00635     if ( (rv = all_mappable_surfs(opts, p)) != 0 )
00636         RETURN(rv);
00637 
00638     if ( opts->debug > 1 )
00639         fprintf(stderr,"-- timing: Surf (all_map)    : time = %.3f\n",
00640                 COX_clock_time());
00641 
00642     if ( (rv = verify_surf_t(opts, p)) != 0)
00643         RETURN(rv);
00644 
00645     if ( opts->debug > 1 )
00646         fprintf(stderr,"-- timing: Surf (verify surf): time = %.3f\n",
00647                 COX_clock_time());
00648 
00649     if ( (rv = get_surf_measures(opts, p)) != 0)
00650         RETURN(rv);
00651 
00652     if ( opts->debug > 1 )
00653         fprintf(stderr,"-- timing: Surf (get measure): time = %.3f\n",
00654                 COX_clock_time());
00655 
00656     RETURN(0);
00657 }

00665 {
00666     int * fcodes;
00667     int   c, debug;
00668     int   geta, getb;
00669 
00670 ENTRY("get_surf_measures");
00671 
00672     geta = getb = 0;
00673     debug = opts->debug > 2;
00674 
00675     if ( opts->info & ST_INFO_VOL )
00676     {
00677         geta = getb = 1;
00678     }
00679     else                /* maybe there are functions that need them */
00680     {
00681         fcodes = p->F->codes;
00682         for ( c = 0; c < p->F->nused; c++ )
00683         {
00684             if ( fcodes[c] == E_SM_N_AREA_A )
00685                 geta = 1;
00686             else if ( fcodes[c] == E_SM_N_AREA_B )
00687                 getb = 1;
00688             else if ( fcodes[c] == E_SM_N_AVEAREA_A )
00689                 geta = 1;
00690             else if ( fcodes[c] == E_SM_N_AVEAREA_B )
00691                 getb = 1;
00692             else if ( fcodes[c] == E_SM_NODE_VOL )
00693             {
00694                 geta = getb = 1;
00695             }
00696         }
00697     }
00698 
00699     if ( geta ) 
00700     {
00701         if ( !SUMA_SurfaceMetrics_eng(p->S.slist[0], "PolyArea", NULL, debug,
00702                                       SUMAg_CF->DsetList) )
00703         {
00704             fprintf(stderr,"** gsf: surface metrics A failure\n");
00705             RETURN(-1);
00706         }
00707 
00708         if ( compute_node_areas(opts, p, 0) != 0 )      /* index 0 */
00709             RETURN(-1);
00710     }
00711 
00712     if ( getb )
00713     {
00714         if ( p->S.nsurf < 2 )
00715         {
00716             fprintf(stderr,"** gsf: functions requre 2 surfaces, failing...\n");
00717             RETURN(-1);
00718         }
00719 
00720         if ( !SUMA_SurfaceMetrics_eng(p->S.slist[1], "PolyArea", NULL, debug,
00721                                       SUMAg_CF->DsetList) )
00722         {
00723             fprintf(stderr,"** gsf: surface metrics B failure\n");
00724             RETURN(-1);
00725         }
00726 
00727         if ( compute_node_areas(opts, p, 1) != 0 )  /* index 1 */
00728             RETURN(-1);
00729     }
00730 
00731     if( geta && getb )
00732     {
00733         if ( surf_triangle_match(opts, p) != 0 )
00734             RETURN(-1);
00735         if ( compute_face_vols(opts, p) != 0 )
00736             RETURN(-1);
00737         if ( compute_node_vols(opts, p) != 0 )
00738             RETURN(-1);
00739     }
00740 
00741     RETURN(0);
00742 }

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     /* clear out the options structure, pointers get explicit NULL */
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;                   /* init to something invalid */
01745     opts->dvox  = -1;                   /* init to something invalid */
01746 
01747     for ( ac = 1; ac < argc; ac++ )
01748     {
01749         /* alphabetical... */
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) )  /* mapping function */
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];     /* store user string for now */
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     /* invalid option */
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 }

01336 {
01337     int c;
01338 
01339 ENTRY("init_opts_t");
01340 
01341     memset(opts, 0, sizeof(opts_t));
01342 
01343     /* init the function list func_t struct */
01344     opts->F.names  = NULL;
01345     opts->F.codes  = NULL;
01346     opts->F.nalloc = 0;
01347     opts->F.nused  = 0;
01348 
01349     /* just to try this out, init the info list */
01350     if ( add_to_flist(&opts->F, g_sm_names[E_SM_NODES]) != 0 )
01351     {
01352         fprintf(stderr,"** failed to init func_t list with '%s'\n",
01353                 g_sm_names[E_SM_NODES]);
01354         RETURN(-1);
01355     }
01356 
01357     opts->spec_file     = NULL;
01358     opts->sv_file       = NULL;
01359     opts->out_1D_file   = NULL;
01360     opts->cmask_cmd     = NULL;
01361     opts->nodes_1D_file = NULL;
01362 
01363     for ( c = 0; c < ST_MAX_SURFS; c++ )
01364         opts->surf_names[c] = NULL;
01365 
01366     opts->dnode         = -1;           /* init to something invalid */
01367 
01368     RETURN(0);
01369 }

02470 {
02471     static int   total_size = 0;
02472     double     * newm, * rowp, * cur, rv;
02473     int          c, row, col, w2, sign;
02474 
02475     if ( width < 2 )
02476         return *data;
02477 
02478     if ( width == 2 )
02479         return( data[0]*data[3] - data[1]*data[2] );
02480 
02481     /* otherwise, we need more space */
02482     w2 = width - 1;
02483     total_size += w2*w2;
02484     newm = malloc(w2*w2 * sizeof(double));
02485     if ( !newm )
02486     {
02487         fprintf(stderr,"** failed to allocate %d doubles for mat (%d ttl)\n",
02488                 w2*w2, total_size);
02489         return 0.0;
02490     }
02491 
02492     rv   = 0.0;
02493     sign = -1;
02494     for ( c = 0; c < width; c++ )
02495     {
02496         sign = -sign;                           /* catch the continue */
02497 
02498         if ( data[c] == 0.0 )                   /* skip any zeros */
02499             continue;
02500 
02501         /* start by copying new data */
02502         for ( row = 1; row < width; row++ )     /* always skip row 0 */
02503         {
02504             rowp = data + row     * width;      /* point to source column */
02505             cur  = newm + (row-1) * w2;         /* init dest pointer */
02506 
02507             for ( col = 0; col < width; col++ )
02508             {
02509                 if ( col == c )
02510                     continue;
02511 
02512                 *cur++ = rowp[col];
02513             }
02514         }
02515 
02516         rv += sign * data[c] * lazy_det( w2, newm );
02517     }
02518 
02519     free(newm);
02520 
02521     return rv;
02522 }

00428 {
00429     float seg[3];
00430     int   c;
00431 
00432 ENTRY("norm2seg_angle");
00433 
00434     for ( c = 0; c < 3; c++ )
00435         seg[c] = p1->xyz[c] - p0->xyz[c];
00436 
00437     RETURN(vector_angle(seg, norm));
00438 }

00498 {
00499     int c, c2, num2print;
00500 
00501 ENTRY("print_column_headers");
00502 
00503     fputc('#', p->outfp);
00504     for (c = 0; c < p->F->nused; c++ )
00505     {
00506         if ( ( p->F->codes[c] == E_SM_COORD_A ) ||
00507              ( p->F->codes[c] == E_SM_COORD_B ) ||
00508              ( p->F->codes[c] == E_SM_NORM_A  ) ||
00509              ( p->F->codes[c] == E_SM_NORM_B  ) )
00510             num2print = 3;
00511         else
00512             num2print = 1;
00513 
00514         for (c2 = 0; c2 < num2print; c2++)
00515             if ( num2print > 1 )
00516                 fprintf(p->outfp, " %8s[%d]", p->F->names[c],c2);
00517             else
00518                 fprintf(p->outfp, "  %10s", p->F->names[c]);
00519     }
00520     fputc('\n', p->outfp);
00521 
00522     fputc('#', p->outfp);
00523     for (c = 0; c < p->F->nused; c++ )
00524     {
00525         if ( ( p->F->codes[c] == E_SM_COORD_A ) ||
00526              ( p->F->codes[c] == E_SM_COORD_B ) ||
00527              ( p->F->codes[c] == E_SM_NORM_A  ) ||
00528              ( p->F->codes[c] == E_SM_NORM_B  ) )
00529             num2print = 3;
00530         else
00531             num2print = 1;
00532 
00533         for (c2 = 0; c2 < num2print; c2++)
00534             fprintf(p->outfp, "  ----------");
00535     }
00536     fputc('\n', p->outfp);
00537 
00538     RETURN(0);
00539 }

01650 {
01651     MRI_IMAGE * im;
01652     float     * fim;
01653     char      * nfile;
01654     int         index;
01655 
01656 ENTRY("read_nodes_file");
01657 
01658     nfile = opts->nodes_1D_file;                /* for ease of typing */
01659 
01660     if ( !nfile )
01661         RETURN(0);
01662 
01663     if ( (im = mri_read_1D(nfile)) == NULL )
01664     {
01665         fprintf(stderr,"** failed to read 1D nodes file '%s'\n", nfile);
01666         RETURN(-1);
01667     }
01668     
01669     if ( im->nx < 1 )
01670     {
01671         fprintf(stderr,"** node list file '%s' appears empty\n", nfile);
01672         RETURN(-1);
01673     }
01674 
01675     if ( im->ny > 1 )
01676     {
01677         fprintf(stderr,"** please specify only one column of node file '%s'\n"
01678                        "   e.g.  -nodes_1D 'surf_data.1D[0]'\n", nfile);
01679         RETURN(-1);
01680     }
01681 
01682     p->nnodes = im->nx;
01683     p->nodes  = (int *)malloc(im->nx * sizeof(int));
01684     ALLOC_CHECK(p->nodes, "int", im->nx);
01685 
01686     /* now get values */
01687     fim = MRI_FLOAT_PTR(im);
01688     for ( index = 0; index < p->nnodes; index++ )
01689         p->nodes[index] = (int)fim[index];
01690 
01691     if ( opts->debug > 0 )
01692         fprintf(stderr,"++ read %d node indices from '%s'\n", p->nnodes, nfile);
01693 
01694     mri_free(im);               /* we're done with image */
01695 
01696     RETURN(0);
01697 }

01222 {
01223     int rv;
01224 
01225 ENTRY("spec2SUMA");
01226 
01227     /* initialize common fields struct */
01228     SUMAg_CF = SUMA_Create_CommonFields();
01229 
01230     if ( SUMAg_CF == NULL )
01231     {
01232         fprintf( stderr, "** failed SUMA_Create_CommonFields(), exiting...\n" );
01233         RETURN(-1);
01234     }
01235 
01236     /* for SUMA type notifications */
01237     if ( opts->debug > 3 )
01238     {
01239         SUMAg_CF->MemTrace = 1;
01240 
01241         if ( opts->debug > 4 )
01242             SUMAg_CF->InOut_Notify = 1;
01243     }
01244 
01245     SUMAg_DOv = SUMA_Alloc_DisplayObject_Struct(SUMA_MAX_DISPLAYABLE_OBJECTS);
01246 
01247     if ( SUMA_Read_SpecFile( opts->spec_file, spec) == 0 )
01248     {
01249         fprintf( stderr, "** failed SUMA_Read_SpecFile(), exiting...\n" );
01250         RETURN(-1);
01251     }
01252 
01253     if ( opts->debug > 2 )
01254         SUMA_ShowSpecStruct(spec, stderr, 3);
01255 
01256     rv = SUMA_spec_select_surfs(spec,opts->surf_names,ST_MAX_SURFS,opts->debug);
01257     if ( rv < 1 )
01258     {
01259         if ( rv == 0 )
01260             fprintf(stderr,"** no given surfaces found in spec file\n");
01261         RETURN(-1);
01262     }
01263 
01264     if ( opts->debug > 1 )
01265         SUMA_ShowSpecStruct(spec, stderr, opts->debug > 2 ? 3 : 1);
01266 
01267     if ( SUMA_spec_set_map_refs(spec, opts->debug) != 0 )
01268         RETURN(-1);
01269 
01270     /* make sure only group was read from spec file */
01271     if ( spec->N_Groups != 1 )
01272     {
01273         fprintf( stderr,"** error: N_Groups <%d> must be 1 in spec file <%s>\n",
01274                  spec->N_Groups, opts->spec_file );
01275         RETURN(-1);
01276     }
01277 
01278     /* actually load the surface(s) from the spec file */
01279     if (SUMA_LoadSpec_eng(spec, SUMAg_DOv, &SUMAg_N_DOv, opts->sv_file,
01280                           opts->debug>3, SUMAg_CF->DsetList) == 0)
01281     {
01282         fprintf( stderr, "** error: failed SUMA_LoadSpec(), exiting...\n" );
01283         RETURN(-1);
01284     }
01285 
01286     if ( opts->debug > 0 )
01287         fputs( "++ surfaces loaded.\n", stderr );
01288 
01289     RETURN(0);
01290 }

00750 {
00751     int   face, maxf;
00752     int * f0, * f1;
00753 
00754 ENTRY("surf_triangle_match");
00755 
00756     f0 = p->S.slist[0]->FaceSetList;
00757     f1 = p->S.slist[1]->FaceSetList;
00758 
00759     maxf = p->S.slist[0]->N_FaceSet * 3;
00760 
00761     if ( !f0 || !f1 || maxf <= 0 )
00762     {
00763         fprintf(stderr,"** stm: bad face data (%p, %p, %d)\n", f0, f1, maxf);
00764         RETURN(-1);
00765     }
00766 
00767     for ( face = 0; face < maxf; face++ )
00768     {
00769         if (*f0 != *f1 )
00770         {
00771             fprintf(stderr,"** face diff @ f3 = %d, *f0,*f1 = (%d,%d)\n",
00772                     face, *f0, *f1);
00773             RETURN(-1);
00774         }
00775         f0++;  f1++;
00776     }
00777 
00778     if ( opts->debug > 1 )
00779         fprintf(stderr,"-- faces okay, woohoo!\n");
00780 
00781     RETURN(0);
00782 }

02266 {
02267     double v0[3], v1[3], v2[3], xp[3];
02268 
02269     v0[0] = p1[0] - p0[0];
02270     v0[1] = p1[1] - p0[1];
02271     v0[2] = p1[2] - p0[2];
02272 
02273     v1[0] = p2[0] - p0[0];
02274     v1[1] = p2[1] - p0[1];
02275     v1[2] = p2[2] - p0[2];
02276 
02277     v2[0] = p3[0] - p0[0];
02278     v2[1] = p3[1] - p0[1];
02279     v2[2] = p3[2] - p0[2];
02280 
02281     cross_product(xp, v1, v2);
02282 
02283     return(fabs(dot_product(v0,xp))/6.0);
02284 }

00836 {
00837     if ( level == USE_SHORT )
00838     {
00839         printf( "usage: %s -help\n"
00840                 "usage: %s [options] file1 file2 ...\n",
00841                 prog, prog);
00842         return 0;
00843     }
00844     else if ( level == USE_VERSION )
00845     {
00846         printf( "%s, version %s, compiled: %s\n",
00847                 prog, VERSION, __DATE__ );
00848         return 0;
00849     }
00850     else if ( level == USE_HISTORY )
00851     {
00852         fputs( g_history, stdout );
00853         return 0;
00854     }
00855     else if ( level == USE_GE )
00856     {
00857         help_ge_structs( prog );
00858         return 0;
00859     }
00860     else if ( level != USE_LONG )
00861     {
00862         fprintf( stderr, "failure: bad usage level <%d>\n", level );
00863         return -1;
00864     }
00865 
00866     /* so we are in a USE_LONG case */
00867     help_full( prog );
00868 
00869     return 0;
00870 }

01715 {
01716     int * fcodes, c, errs;
01717 
01718 ENTRY("validate_option_lists");
01719 
01720     errs = 0;
01721 
01722     if ( (opts->info & ST_INFO_THICK) && p->S.nsurf < 2 )
01723     {
01724         fprintf(stderr,"** -info_thick option requiers 2 surfaces\n");
01725         errs++;
01726     }
01727 
01728     if ( (opts->info & ST_INFO_VOL) && p->S.nsurf < 2 )
01729     {
01730         fprintf(stderr,"** -info_vol option requiers 2 surfaces\n");
01731         errs++;
01732     }
01733 
01734     /* verify all the functions in our list */
01735     fcodes = p->F->codes;       /* just for less typing */
01736     for ( c = 0; c < p->F->nused; c++ )
01737     {
01738         switch (fcodes[c])
01739         {
01740             default:
01741                 fprintf(stderr, "** vol: invalid function code %d\n",fcodes[c]);
01742                 errs++;
01743                 break;
01744 
01745             case E_SM_COORD_A:
01746             case E_SM_N_AREA_A:
01747             case E_SM_N_AVEAREA_A:
01748             case E_SM_NODES:
01749             case E_SM_NORM_A:
01750                 break;
01751 
01752             case E_SM_ANG_NS_A:
01753                 if ( !p->S.slist[0]->NodeNormList )
01754                 {
01755                     fprintf(stderr,"** missing node normals for func '%s'\n",
01756                             g_sm_names[fcodes[c]]);
01757                     errs++;
01758                 }
01759                 break;
01760 
01761             case E_SM_ANG_NORMS:
01762             case E_SM_ANG_NS_B:
01763             case E_SM_NORM_B:
01764                 if ( !p->S.slist[0]->NodeNormList ||
01765                      !p->S.slist[1]->NodeNormList )
01766                 {
01767                     fprintf(stderr,"** missing node normals for func '%s'\n",
01768                             g_sm_names[fcodes[c]]);
01769                     errs++;
01770                 }
01771 
01772                 SM_2SURF_TEST(c);
01773                 break;
01774                 
01775             case E_SM_COORD_B:
01776             case E_SM_N_AREA_B:
01777             case E_SM_N_AVEAREA_B:
01778             case E_SM_NODE_VOL:
01779             case E_SM_THICK:
01780 
01781                 SM_2SURF_TEST(c);
01782                 break;
01783 
01784             case E_SM_NTRI:
01785                 if ( ! p->S.slist[0]->MF->N_Memb )
01786                 {
01787                     fprintf(stderr,"** missing Face Member list, func '%s'\n",
01788                             g_sm_names[fcodes[c]]);
01789                     errs++;
01790                 }
01791 
01792                 break;
01793         }
01794 
01795         if ( opts->debug > 1 && errs == 0 )
01796             fprintf(stderr,"++ have valid function name '%s'\n",
01797                     g_sm_names[fcodes[c]]);
01798     }
01799 
01800     if ( opts->debug > 0 && errs > 0 )
01801         fprintf(stderr,"** validate_option_lists: %d errors found\n", errs);
01802 
01803     if ( errs > 0 )
01804         RETURN(-1);
01805 
01806     if ( opts->debug > 0 )
01807         fprintf(stderr,"-- found %d output functions\n", p->F->nused);
01808 
01809     RETURN(0);
01810 }

02042 {
02043 ENTRY("validate_options");
02044 
02045     /* clear param struct - pointer get explicit NULL */
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 }

00468 {
00469     double mag0, mag1, dot, ratio, angle;
00470 
00471 ENTRY("vector_angle");
00472 
00473     mag0 = fvec_magnitude(v0, 3);
00474     mag1 = fvec_magnitude(v1, 3);
00475     dot  = v0[0]*v1[0] + v0[1]*v1[1] + v0[2]*v1[2];
00476 
00477     if ( mag0 == 0.0 || mag1 == 0.0 )
00478         RETURN(0.0);
00479 
00480     ratio = dot / (mag0 * mag1);
00481     RANGE(-1.0, ratio, 1.0);
00482 
00483     angle = acos(ratio);
00484 
00485     /* keep angle in [0,PI/2] */
00486     if ( 2 * angle > ST_PI )
00487         angle = PI - angle;
00488 
00489     RETURN(angle);
00490 }

00548 {
00549     int c;
00550 
00551 ENTRY("verify_surf_t");
00552 
00553     if ( p->S.nsurf < 1 )
00554     {
00555         fprintf(stderr,"** no surfaces found\n");
00556         if ( opts->debug > 0 )
00557             disp_surf_t("-- empty? ", &p->S);
00558         RETURN(-1);
00559     }
00560 
00561     if ( p->S.nsurf > 1 )
00562     {
00563         for ( c = 1; c < p->S.nsurf; c++ )
00564             if ( p->S.slist[c]->N_Node != p->S.nnodes )
00565             {
00566                 fprintf(stderr,"** surface %d has %d nodes (should be %d)\n",
00567                         c, p->S.slist[c]->N_Node, p->S.nnodes);
00568                 RETURN(-1);
00569             }
00570     }
00571 
00572     /* verify nodes and cmask */
00573     if ( ! p->nodes )                   /* if empty, create the trivial list */
00574     {
00575         p->nnodes = p->S.nnodes;
00576         p->nodes  = (int *)malloc(p->nnodes * sizeof(int));
00577         ALLOC_CHECK(p->nodes, "int", p->S.nnodes);
00578 
00579         /* now fill the list with trivial indices */
00580         for ( c = 0; c < p->nnodes; c++ )
00581             p->nodes[c] = c;
00582     }
00583     else                                /* verify that the indices are valid */
00584     {
00585         for ( c = 0; c < p->nnodes; c++ )
00586             if ( (p->nodes[c] < 0) || (p->nodes[c] >= p->S.nnodes) )
00587             {
00588                 fprintf(stderr,"** error: node list index %d (value = %d):"
00589                                "          outside valid range [0,%d]\n",
00590                         c, p->nodes[c], p->S.nnodes-1);
00591                 RETURN(-1);
00592             }
00593     }
00594 
00595     /* in any case, if there is a mask, the length should match nnodes */
00596     if ( p->cmask && (p->ncmask != p->nnodes) )
00597     {
00598         fprintf(stderr,"** cmask and node list lengths differ (%d, %d)\n",
00599                 p->ncmask, p->nnodes);
00600         RETURN(-1);
00601     }
00602 
00603     if ( opts->debug > 1 )
00604     {
00605         disp_param_t("-- surf params verified: ", p);
00606         disp_surf_t ("-- surf params verified: ", &p->S);
00607     }
00608 
00609     if ( validate_option_lists(opts, p) != 0 )
00610         RETURN(-1);
00611 
00612     RETURN(0);
00613 }

00197 {
00198     SUMA_SurfaceObject * so0, * so1;
00199 
00200     THD_fvec3   p0, p1;
00201     double      tarea0, tarea1, tvolume;
00202     double      dist, min_dist, max_dist, tdist;
00203     double      atn = 0.0, atna = 0.0, atnb = 0.0;  /* angle totals */
00204     float       fvn, fva, fvb;
00205     float     * fp0, * fp1;
00206     float     * norms0, * norms1;
00207     float       ave_dist, fval;
00208     int         c, fnum, node, nindex;
00209     int       * fcodes;
00210     int         skipped = 0;
00211 
00212 ENTRY("write_output");
00213 
00214     so0 = p->S.slist[0];
00215     so1 = p->S.slist[1];
00216 
00217     /* just a simple case for now */
00218 
00219     print_column_headers(opts, p);
00220 
00221     /* initialize some pointers */
00222     fp0    = so0->NodeList;
00223     norms0 = so0->NodeNormList;
00224 
00225     /* if we have only one surface, just init fp1 to that */
00226     if ( p->S.nsurf > 1 )
00227     {
00228         fp1    = so1->NodeList;
00229         norms1 = so1->NodeNormList;
00230     }
00231     else
00232     {
00233         fp1    = fp0;
00234         norms1 = norms0;
00235     }
00236 
00237     fcodes = p->F->codes;       /* for convenience */
00238 
00239     tdist   = 0.0;              /* for total distance over nodes */
00240     tarea0  = 0.0;              /* for total area computation   */
00241     tarea1  = 0.0;              /* for total area computation  */
00242     tvolume = 0.0;              /* for total volume           */
00243 
00244     min_dist = 9999.0;
00245     max_dist = 0.0;
00246     for (nindex = 0; nindex < p->nnodes; nindex++)
00247     {
00248         if ( p->cmask && !p->cmask[nindex] )
00249         {
00250             skipped++;
00251             continue;
00252         }
00253 
00254         node = p->nodes[nindex];
00255 
00256         memcpy(p0.xyz, fp0+3*node, 3*sizeof(float));
00257         memcpy(p1.xyz, fp1+3*node, 3*sizeof(float));
00258 
00259         dist = dist_fn(3, p0.xyz, p1.xyz);
00260         if ( dist < min_dist ) min_dist = dist;
00261         if ( dist > max_dist ) max_dist = dist;
00262         tdist += dist;
00263 
00264         /* keep track of the total areas and volume */
00265         if ( p->S.narea[0] )
00266             tarea0  += p->S.narea[0][node];
00267 
00268         if ( p->S.narea[1] )
00269             tarea1  += p->S.narea[1][node];
00270 
00271         if ( p->S.nvol )
00272             tvolume += p->S.nvol[node];
00273 
00274         fputc(' ', p->outfp);
00275 
00276         for (fnum = 0; fnum < p->F->nused; fnum++)
00277         {
00278             switch(fcodes[fnum])
00279             {
00280                 default:
00281                     fprintf(stderr,"** bad output Info code %d\n",fcodes[fnum]);
00282                     break;
00283 
00284                 case E_SM_ANG_NORMS:
00285                     fvn = vector_angle(norms0 + 3*node, norms1 + 3*node);
00286                     fprintf(p->outfp,"  %10s", MV_format_fval(fvn));
00287                     atn += fvn;
00288                     break;
00289 
00290                 case E_SM_ANG_NS_A:
00291                     fva = norm2seg_angle(&p0, &p1, norms0 + 3*node);
00292                     fprintf(p->outfp,"  %10s", MV_format_fval(fva));
00293                     atna += fva;
00294                     break;
00295 
00296                 case E_SM_ANG_NS_B:
00297                     fvb = norm2seg_angle(&p0, &p1, norms1 + 3*node);
00298                     fprintf(p->outfp,"  %10s", MV_format_fval(fvb));
00299                     atnb += fvb;
00300                     break;
00301 
00302                 case E_SM_COORD_A:
00303                     for (c = 0; c < 3; c++)
00304                         fprintf(p->outfp,"  %10s",
00305                                 MV_format_fval(fp0[3*node+c]));
00306                     break;
00307 
00308                 case E_SM_COORD_B:
00309                     for (c = 0; c < 3; c++)
00310                         fprintf(p->outfp,"  %10s",
00311                                 MV_format_fval(fp1[3*node+c]));
00312                     break;
00313 
00314                 case E_SM_N_AREA_A:
00315                     fprintf(p->outfp,"  %10s",
00316                             MV_format_fval(p->S.narea[0][node]));
00317                     break;
00318 
00319                 case E_SM_N_AREA_B:
00320                     fprintf(p->outfp,"  %10s",
00321                             MV_format_fval(p->S.narea[1][node]));
00322                     break;
00323 
00324                 case E_SM_N_AVEAREA_A:
00325                     fval = 3 * p->S.narea[0][node] / so0->MF->N_Memb[node];
00326                     fprintf(p->outfp,"  %10s", MV_format_fval(fval));
00327                     break;
00328 
00329                 case E_SM_N_AVEAREA_B:
00330                     fval = 3 * p->S.narea[1][node] / so1->MF->N_Memb[node];
00331                     fprintf(p->outfp,"  %10s", MV_format_fval(fval));
00332                     break;
00333 
00334                 case E_SM_NTRI:
00335                     fprintf(p->outfp,"  %10d", so0->MF->N_Memb[node]);
00336                     break;
00337 
00338                 case E_SM_NODES:
00339                     fprintf(p->outfp,"  %10d", node);
00340                     break;
00341 
00342                 case E_SM_NODE_VOL:
00343                     fprintf(p->outfp,"  %10s", MV_format_fval(p->S.nvol[node]));
00344                     break;
00345 
00346                 case E_SM_NORM_A:
00347                     if ( norms0 )
00348                         for (c = 0; c < 3; c++)
00349                             fprintf(p->outfp,"  %10s",
00350                                     MV_format_fval(norms0[3*node+c]));
00351                     break;
00352 
00353                 case E_SM_NORM_B:
00354                     if ( norms1 )
00355                         for (c = 0; c < 3; c++)
00356                             fprintf(p->outfp,"  %10s",
00357                                     MV_format_fval(norms1[3*node+c]));
00358                     break;
00359 
00360                 case E_SM_THICK:
00361                     fprintf(p->outfp,"  %10s", MV_format_fval(dist));
00362                     break;
00363             }
00364         }
00365 
00366         fputc('\n', p->outfp);
00367 
00368         if ( node == opts->dnode && opts->debug > 0 )
00369         {
00370             fprintf(stderr,"-- dnode %d:\n", node);
00371             disp_f3_point("-- p0    = ", fp0+3*node);
00372             disp_f3_point("-- p1    = ", fp1+3*node);
00373 
00374             if ( atn > 0.0 || atna > 0.0 || atnb > 0.0 )
00375             {
00376                 disp_f3_point("-- normA = ", norms0 + 3*node);
00377                 disp_f3_point("-- normB = ", norms1 + 3*node);
00378             }
00379         }
00380     }
00381 
00382     ave_dist = tdist/p->nnodes;
00383 
00384     if ( opts->info )
00385         printf("----------------------------------------------------------\n");
00386 
00387     if ( opts->info & ST_INFO_AREA )
00388     {
00389         if ( p->S.narea[0] )
00390             printf("-- total area 0 = %.1f\n", tarea0);
00391 
00392         if ( p->S.nsurf > 1 )
00393             printf("-- total area 1 = %.1f\n", tarea1);
00394     }
00395 
00396     if ( (opts->info & ST_INFO_NORMS) &&
00397          (atn > 0.0 || atna > 0.0 || atnb > 0.0) )
00398     {
00399         printf( "-- ave. angles to normals: (nA_nB, sA, sB) = "
00400                 "(%.4f, %.4f, %.4f)\n",
00401                 atn/p->nnodes, atna/p->nnodes, atnb/p->nnodes);
00402     }
00403 
00404     if ( (opts->info & ST_INFO_THICK) && (p->S.nsurf > 1) )
00405         printf( "-- thickness: (min, ave, max) = (%.5f, %.5f, %.5f)\n",
00406                 min_dist, ave_dist, max_dist);
00407 
00408     if ( (opts->info & ST_INFO_VOL) && p->S.nvol )
00409     {
00410         printf("-- total volume = %.1f\n", tvolume);
00411         if ( opts->debug > 1 )
00412             printf("-- ave dist * (area0, ave, area1) = (%.1f, %.1f, %.1f)\n",
00413                 ave_dist*tarea0, ave_dist*(tarea0+tarea1)/2, ave_dist*tarea1);
00414     }
00415 
00416     if ( p->cmask )
00417         printf("-- from cmask, nodes skipped = %d\n", skipped);
00418 
00419     RETURN(0);
00420 }