00001
00002 #define VERSION "version 1.11 (October 6, 2004)"
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045 static char g_history[] =
00046 "----------------------------------------------------------------------\n"
00047 "history :\n"
00048 "\n"
00049 "1.0 December 01, 2003 [rickr]\n"
00050 " - initial release\n"
00051 "\n"
00052 "1.1 December 02, 2003 [rickr]\n"
00053 " - fixed stupid macro error, grrrr...\n"
00054 "\n"
00055 "1.2 December 03, 2003 [rickr]\n"
00056 " - added '-cmask' and '-nodes_1D' options\n"
00057 "\n"
00058 "1.3 December 11, 2003 [rickr]\n"
00059 " - added required program argument(s): '-surf_A' (and 'B' for 2 surfs)\n"
00060 " o see '-help' for a description and examples\n"
00061 " - added SUMA_spec_select_surfs() and SUMA_spec_set_map_refs() so that:\n"
00062 " o only requested surfaces are loaded\n"
00063 " o spec files no longer need 'SAME' as MappingRef\n"
00064 " - fixed loss of default node indices (from nodes_1D change)\n"
00065 " - added '-hist' option\n"
00066 " - display angle averages only if at least 1 total is computed\n"
00067 "\n"
00068 "1.4 January 22, 2004 [rickr]\n"
00069 " - fixed error with '-nodes_1D' indexing (fp0 in write_output())\n"
00070 " (for output of node coordinates)\n"
00071 " - added '-sv' option to examples\n"
00072 " - reversed history list (most recent last) for '-hist' option\n"
00073 "\n"
00074 "1.5 January 23, 2004 [rickr]\n"
00075 " - SUMA_isINHmappable() is deprecated, check with AnatCorrect field\n"
00076 "\n"
00077 "1.6 February 11, 2004 [rickr]\n"
00078 " - add a little debug help for !AnatCorrect case\n"
00079 "\n"
00080 "1.7 February 23, 2004 [rickr]\n"
00081 " - added functions 'n_avearea_A', 'n_avearea_B', 'n_ntri'\n"
00082 "\n"
00083 "1.8 March 26, 2004 [ziad]\n"
00084 " - DsetList added to SUMA_LoadSpec_eng() and SUMA_SurfaceMetrics_eng()\n"
00085 "\n"
00086 "1.9 July 29, 2004 [rickr]\n"
00087 " - Remove check for anat correct.\n"
00088 "\n"
00089 "1.10 August 11, 2004 [rickr]\n"
00090 " - Add comment about volume being too large.\n"
00091 "\n"
00092 "1.11 October 12, 2004 [rickr]\n"
00093 " - More volume comments: suggest 'SurfPatch -vol'\n"
00094 "----------------------------------------------------------------------\n";
00095
00096
00097
00098
00099
00100
00101 #include "mrilib.h"
00102 #include "SUMA_suma.h"
00103 #include "SUMA_SurfMeasures.h"
00104
00105 extern void machdep( void );
00106
00107
00108
00109 SUMA_SurfaceViewer * SUMAg_SVv = NULL;
00110 int SUMAg_N_SVv = 0;
00111 SUMA_DO * SUMAg_DOv = NULL;
00112 int SUMAg_N_DOv = 0;
00113 SUMA_CommonFields * SUMAg_CF = NULL;
00114
00115
00116 char * g_sm_names[] = { "none", "ang_norms", "ang_ns_A", "ang_ns_B",
00117 "coord_A", "coord_B", "n_area_A", "n_area_B",
00118 "n_avearea_A", "n_avearea_B", "n_ntri",
00119 "node_vol", "nodes", "norm_A", "norm_B", "thick" };
00120
00121 char * g_sm_desc[] = { "invalid function",
00122 "angular difference between normals",
00123 "angular diff between segment and first norm",
00124 "angular diff between segment and second norm",
00125 "xyz coordinates of node on first surface",
00126 "xyz coordinates of node on second surface",
00127 "associated node area on first surface",
00128 "associated node area on second surface",
00129 "for each node, average area of triangles (surf A)",
00130 "for each node, average area of triangles (surf B)",
00131 "for each node, number of associated triangles",
00132 "associated node volume between surfaces",
00133 "node number",
00134 "vector of normal at node on first surface",
00135 "vector of normal at node on second surface",
00136 "distance between surfaces along segment" };
00137
00138
00139
00140 int main ( int argc, char * argv[] )
00141 {
00142 param_t p;
00143 opts_t opts;
00144 int rv;
00145
00146
00147 (void) COX_clock_time();
00148
00149 mainENTRY("SurfMeasures main");
00150 machdep();
00151 AFNI_logger(PROG_NAME,argc,argv);
00152
00153 if ( (rv = init_options(&opts, argc, argv)) != 0 )
00154 return rv;
00155
00156 if ( opts.debug > 1 )
00157 fprintf(stderr,"-- timing: init opts : time = %.3f\n",
00158 COX_clock_time());
00159
00160 if ( (rv = validate_options(&opts, &p)) != 0 )
00161 return rv;
00162
00163 if ( opts.debug > 1 )
00164 fprintf(stderr,"-- timing: validate opts : time = %.3f\n",
00165 COX_clock_time());
00166
00167 if ( (rv = get_surf_data(&opts, &p)) != 0 )
00168 return rv;
00169
00170 if ( opts.debug > 1 )
00171 fprintf(stderr,"-- timing: get surf data : time = %.3f\n",
00172 COX_clock_time());
00173
00174 if ( (rv = write_output(&opts, &p)) != 0 )
00175 return rv;
00176
00177 if ( opts.debug > 1 )
00178 fprintf(stderr,"-- timing: write outfile : time = %.3f\n",
00179 COX_clock_time());
00180
00181 final_cleanup(&opts, &p);
00182
00183 if ( opts.debug > 1 )
00184 fprintf(stderr,"-- timing: final clean up : time = %.3f\n",
00185 COX_clock_time());
00186
00187 return 0;
00188 }
00189
00190
00191
00192
00193
00194
00195
00196 int write_output( opts_t * opts, param_t * p )
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;
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
00218
00219 print_column_headers(opts, p);
00220
00221
00222 fp0 = so0->NodeList;
00223 norms0 = so0->NodeNormList;
00224
00225
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;
00238
00239 tdist = 0.0;
00240 tarea0 = 0.0;
00241 tarea1 = 0.0;
00242 tvolume = 0.0;
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
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 }
00421
00422
00423
00424
00425
00426
00427 float norm2seg_angle( THD_fvec3 * p0, THD_fvec3 * p1, float * norm )
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 }
00439
00440
00441
00442
00443
00444
00445 float fvec_magnitude( float * v, int length )
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 }
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467 float vector_angle( float * v0, float * v1 )
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
00486 if ( 2 * angle > ST_PI )
00487 angle = PI - angle;
00488
00489 RETURN(angle);
00490 }
00491
00492
00493
00494
00495
00496
00497 int print_column_headers( opts_t * opts, param_t * p )
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 }
00540
00541
00542
00543
00544
00545
00546
00547 int verify_surf_t( opts_t * opts, param_t * p )
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
00573 if ( ! p->nodes )
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
00580 for ( c = 0; c < p->nnodes; c++ )
00581 p->nodes[c] = c;
00582 }
00583 else
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
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 }
00614
00615
00616
00617
00618
00619
00620
00621 int get_surf_data( opts_t * opts, param_t * p )
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 }
00658
00659
00660
00661
00662
00663
00664 int get_surf_measures( opts_t * opts, param_t * p )
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
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 )
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 )
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 }
00743
00744
00745
00746
00747
00748
00749 int surf_triangle_match( opts_t * opts, param_t * p )
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 }
00783
00784
00785 #if 0
00786
00787
00788
00789
00790
00791 int test_planar_sides( opts_t * opts, param_t * p )
00792 {
00793 SUMA_SurfaceObject * so;
00794 float * nodes0, * nodes1;
00795 float p1, p2, p3, max, mold;
00796 int a, b, c, mcount;
00797 int face, nfaces;
00798 int * flist;
00799
00800 ENTRY("test_planar_sides");
00801
00802 flist = p->S.slist[0]->FaceSetList;
00803 nfaces = p->S.slist[1]->N_FaceSet;
00804
00805 nodes0 = p->S.slist[0]->NodeList;
00806 nodes1 = p->S.slist[1]->NodeList;
00807
00808 if ( !flist || !nodes0 || !nodes1 || nfaces <= 0 )
00809 {
00810 fprintf(stderr,"** tps: bad face or node data (%p, %p, %p, %d)\n",
00811 flist, nodes0, nodes1, nfaces);
00812 RETURN(-1);
00813 }
00814
00815 max = 0.0;
00816 mold = 1.0;
00817 mcount = 0;
00818 for ( face = 0; face < nfaces; face++ )
00819 {
00820 a = flist[0] * 3;
00821 b = flist[1] * 3;
00822 c = flist[2] * 3;
00823
00824 p1 = eval_planar_points(nodes0+a, nodes0+b, nodes1+a, nodes1+b);
00825 p2 = eval_planar_points(nodes0+b, nodes0+c, nodes1+b, nodes1+c);
00826 p3 = eval_planar_points(nodes0+c, nodes0+a, nodes1+c, nodes1+a);
00827 if (p1 > max) max = p1;
00828 if (p2 > max) max = p2;
00829 if (p3 > max) max = p3;
00830
00831 if ( max > mold )
00832 {
00833 mcount++;
00834 if ( opts->debug > 2 )
00835 {
00836 fprintf(stderr,"** new max (%d) for face %d, m = %f, m0 = %f\n",
00837 mcount, face, max, mold );
00838 disp_f3_point(" a0 = ", nodes0+a);
00839 disp_f3_point(" b0 = ", nodes0+b);
00840 disp_f3_point(" c0 = ", nodes0+c);
00841 disp_f3_point(" a1 = ", nodes1+a);
00842 disp_f3_point(" b1 = ", nodes1+b);
00843 disp_f3_point(" c1 = ", nodes1+c);
00844 }
00845 mold = max*1.5;
00846 }
00847
00848 flist += 3;
00849 }
00850
00851 if ( max > 1.0 )
00852 RETURN(1);
00853
00854 RETURN(0);
00855 }
00856
00857
00858
00859
00860
00861
00862
00863
00864 float eval_planar_points( float * a, float * b, float * c, float * d )
00865 {
00866 ENTRY("planar_points");
00867
00868 double d0[3], d1[3], d2[3];
00869 double xp[3];
00870 double result;
00871
00872 d0[0] = c[0] - a[0];
00873 d0[1] = c[1] - a[1];
00874 d0[2] = c[2] - a[2];
00875
00876 d1[0] = b[0] - a[0];
00877 d1[1] = b[1] - a[1];
00878 d1[2] = b[2] - a[2];
00879
00880 d2[0] = d[0] - c[0];
00881 d2[1] = d[1] - c[1];
00882 d2[2] = d[2] - c[2];
00883
00884 cross_product(xp, d1, d2);
00885
00886 result = d0[0]*xp[0] + d0[1]*xp[1] + d0[2]*xp[2];
00887
00888 RETURN(fabs(result));
00889 }
00890
00891 #endif
00892
00893
00894
00895
00896
00897
00898 int cross_product( double * res, double * u, double * v )
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 }
00906
00907
00908
00909
00910
00911
00912 double dot_product( double * u, double * v )
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 }
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931 int compute_node_areas( opts_t * opts, param_t * p, int sindex )
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];
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 }
00986
00987
00988
00989
00990
00991
00992
00993
00994
00995 int compute_node_vols( opts_t * opts, param_t * p )
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];
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 }
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053 int compute_face_vols( opts_t * opts, param_t * p )
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;
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;
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 }
01146
01147
01148
01149
01150
01151
01152
01153 int all_mappable_surfs( opts_t * opts, param_t * p )
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;
01165 p->S.nsurf = 0;
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
01175
01176 #if 0
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 }
01214
01215
01216
01217
01218
01219
01220
01221 int spec2SUMA( SUMA_SurfSpecFile * spec, opts_t * opts )
01222 {
01223 int rv;
01224
01225 ENTRY("spec2SUMA");
01226
01227
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
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
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
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 }
01291
01292
01293
01294
01295
01296
01297 int add_to_flist( func_t * F, char * fname )
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++;
01326
01327 RETURN(0);
01328 }
01329
01330
01331
01332
01333
01334
01335 int init_opts_t( opts_t * opts )
01336 {
01337 int c;
01338
01339 ENTRY("init_opts_t");
01340
01341 memset(opts, 0, sizeof(opts_t));
01342
01343
01344 opts->F.names = NULL;
01345 opts->F.codes = NULL;
01346 opts->F.nalloc = 0;
01347 opts->F.nused = 0;
01348
01349
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;
01367
01368 RETURN(0);
01369 }
01370
01371
01372
01373
01374
01375 #define CHECK_ARG_COUNT(ac,str) \
01376 do { \
01377 if ((ac+1) >= argc) { \
01378 fputs(str,stderr); \
01379 RETURN(-1); \
01380 } \
01381 } while (0)
01382
01383
01384
01385
01386
01387 int init_options( opts_t * opts, int argc, char * argv[] )
01388 {
01389 int ac, ind;
01390
01391 ENTRY("init_options");
01392
01393 if ( argc < 2 )
01394 RETURN( usage(PROG_NAME, ST_USE_LONG) );
01395
01396
01397 if ( init_opts_t(opts) != 0 )
01398 RETURN(-1);
01399
01400 for ( ac = 1; ac < argc; ac++ )
01401 {
01402 if ( ! strncmp(argv[ac], "-debug", 6) )
01403 {
01404 CHECK_ARG_COUNT(ac,"option usage: -debug LEVEL\n");
01405
01406 opts->debug = atoi(argv[++ac]);
01407 if ( opts->debug < 0 || opts->debug > ST_DEBUG_MAX_LEVEL )
01408 {
01409 fprintf(stderr,"** bad debug level %d, should be in [0,%d]\n",
01410 opts->debug, ST_DEBUG_MAX_LEVEL);
01411 RETURN(-1);
01412 }
01413 }
01414 else if ( ! strncmp(argv[ac], "-cmask", 6) )
01415 {
01416 CHECK_ARG_COUNT(ac,"option usage: -cmask COMMAND\n");
01417 opts->cmask_cmd = argv[++ac];
01418 }
01419 else if ( ! strncmp(argv[ac], "-dnode", 6) )
01420 {
01421 CHECK_ARG_COUNT(ac,"option usage: -dnode NODE_NUM\n");
01422 opts->dnode = atoi(argv[++ac]);
01423 }
01424 else if ( ! strncmp(argv[ac], "-func", 5) )
01425 {
01426 CHECK_ARG_COUNT(ac,"option usage: -func FUNCTION\n");
01427 ++ac;
01428 if ( add_to_flist(&opts->F, argv[ac]) != 0 )
01429 RETURN(-1);
01430 }
01431 else if ( ! strncmp(argv[ac], "-help", 5) )
01432 RETURN(usage(PROG_NAME, ST_USE_LONG));
01433 else if ( ! strncmp(argv[ac], "-hist", 5) )
01434 RETURN(usage(PROG_NAME, ST_USE_HIST));
01435 else if ( ! strncmp(argv[ac], "-info_all",9) )
01436 {
01437 opts->info |= ST_INFO_ALL;
01438 }
01439 else if ( ! strncmp(argv[ac], "-info_area",10) )
01440 {
01441 opts->info |= ST_INFO_AREA;
01442 }
01443 else if ( ! strncmp(argv[ac], "-info_norms",10) )
01444 {
01445 opts->info |= ST_INFO_NORMS;
01446 }
01447 else if ( ! strncmp(argv[ac], "-info_thick",11) )
01448 {
01449 opts->info |= ST_INFO_THICK;
01450 }
01451 else if ( ! strncmp(argv[ac], "-info_vol",9) )
01452 {
01453 opts->info |= ST_INFO_VOL;
01454 }
01455 else if ( ! strncmp(argv[ac], "-nodes_1D", 9) )
01456 {
01457 CHECK_ARG_COUNT(ac,"option usage: -nodes_1D NODE_LIST_FILE\n");
01458 opts->nodes_1D_file = argv[++ac];
01459 }
01460 else if ( ! strncmp(argv[ac], "-out_1D", 7) )
01461 {
01462 CHECK_ARG_COUNT(ac,"option usage: -out_1D OUTPUT_FILE\n");
01463 opts->out_1D_file = argv[++ac];
01464 }
01465 else if ( ! strncmp(argv[ac], "-spec", 5) )
01466 {
01467 CHECK_ARG_COUNT(ac,"option usage: -spec SPEC_FILE\n");
01468 opts->spec_file = argv[++ac];
01469 }
01470 else if ( ! strncmp(argv[ac], "-surf_", 6) )
01471 {
01472 CHECK_ARG_COUNT(ac,"option usage: -surf_X SURF_NAME\n");
01473 ind = argv[ac][6] - 'A';
01474 if ( (ind < 0) || (ind >= ST_MAX_SURFS) )
01475 {
01476 fprintf(stderr,"** -surf_X option '%s' out of range,\n"
01477 " use one of '-surf_A' through '-surf_%c'\n",
01478 argv[ac], 'A'+ST_MAX_SURFS-1);
01479 RETURN(-1);
01480 }
01481 opts->surf_names[ind] = argv[++ac];
01482 }
01483 else if ( ! strncmp(argv[ac], "-sv", 3) )
01484 {
01485 CHECK_ARG_COUNT(ac,"option usage: -sv SURF_VOLUME\n");
01486 opts->sv_file = argv[++ac];
01487 }
01488 else if ( ! strncmp(argv[ac], "-ver",4) )
01489 {
01490 RETURN( usage(PROG_NAME, ST_USE_VERSION) );
01491 }
01492 else
01493 {
01494 fprintf(stderr,"invalid option <%s>\n",argv[ac]);
01495 RETURN( usage(PROG_NAME, ST_USE_SHORT) );
01496 }
01497 }
01498
01499 RETURN(0);
01500 }
01501
01502
01503
01504
01505
01506
01507 int validate_options( opts_t * opts, param_t * p )
01508 {
01509 int errs = 0;
01510
01511 ENTRY("validate_options");
01512
01513 if ( !opts || !p )
01514 {
01515 fprintf(stderr,"** vo: bad params (%p,%p)\n", opts, p);
01516 RETURN(-1);
01517 }
01518
01519 if ( opts->F.nused <= 0 )
01520 {
01521 fprintf(stderr,"** must specify at least one '-func' option\n");
01522 errs++;
01523 }
01524
01525 if ( ! opts->spec_file )
01526 {
01527 fprintf(stderr,"** missing argument: -spec\n");
01528 errs++;
01529 }
01530
01531 if ( ! opts->surf_names[0] )
01532 {
01533 fprintf(stderr,"** missing argument -surf_A\n");
01534 errs++;
01535 }
01536
01537
01538
01539
01540 if ( ! opts->out_1D_file )
01541 {
01542 fprintf(stderr,"** missing argument: -out_1D\n");
01543 errs++;
01544 }
01545 else if ( THD_is_file(opts->out_1D_file) )
01546 {
01547 fprintf(stderr,"** output file already exists: %s\n",opts->out_1D_file);
01548 errs++;
01549 }
01550
01551 if ( errs > 0 )
01552 RETURN(-1);
01553
01554 if ( opts->debug > 1 )
01555 {
01556 disp_opts_t( "-- opts okay: ", opts );
01557 disp_func_t( "-- opts okay: ", &opts->F );
01558 }
01559
01560
01561
01562 memset(p, 0, sizeof(param_t));
01563
01564 p->S.slist = NULL;
01565 p->S.narea[0] = NULL;
01566 p->S.narea[1] = NULL;
01567 p->S.nvol = NULL;
01568 p->S.fvol = NULL;
01569
01570 p->F = &opts->F;
01571
01572 if ( (p->outfp = fopen(opts->out_1D_file, "w")) == NULL )
01573 {
01574 fprintf(stderr,"** cannot open output file '%s'\n",opts->out_1D_file);
01575 RETURN(-1);
01576 }
01577
01578
01579 p->nodes = NULL;
01580 p->nnodes = 0;
01581 p->cmask = NULL;
01582 p->ncmask = 0;
01583
01584 if ( opts->nodes_1D_file )
01585 {
01586 if ( read_nodes_file(opts, p) != 0 )
01587 RETURN(-1);
01588 }
01589
01590 if ( opts->cmask_cmd )
01591 {
01592 if ( get_cmask(opts, p) != 0 )
01593 RETURN(-1);
01594 }
01595
01596 RETURN(0);
01597 }
01598
01599
01600
01601
01602
01603
01604 int get_cmask( opts_t * opts, param_t * p )
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);
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 )
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 }
01643
01644
01645
01646
01647
01648
01649 int read_nodes_file( opts_t * opts, param_t * p )
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;
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
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);
01695
01696 RETURN(0);
01697 }
01698
01699
01700
01701
01702 #define SM_2SURF_TEST(code) \
01703 do{ if (p->S.nsurf<2) { \
01704 fprintf(stderr,"** function %s requires 2 surfaces\n", \
01705 g_sm_names[fcodes[code]]); \
01706 errs++; \
01707 } \
01708 } while (0);
01709
01710
01711
01712
01713
01714 int validate_option_lists( opts_t * opts, param_t * p )
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
01735 fcodes = p->F->codes;
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 }
01811
01812
01813
01814
01815
01816
01817
01818
01819
01820
01821
01822
01823
01824 int usage( char * prog, int use_type )
01825 {
01826 int c;
01827
01828 ENTRY("usage");
01829
01830 if ( use_type == ST_USE_SHORT )
01831 {
01832 fprintf(stderr,"usage: %s [options] -spec SPEC_FILE -func FUNC_NAME\\\n"
01833 " -out_1D OUTFILE\n", prog);
01834 }
01835 else if ( use_type == ST_USE_LONG )
01836 {
01837 printf(
01838 "\n"
01839 "%s - compute measures from the surface dataset(s)\n"
01840 "\n"
01841 " usage: %s [options] -spec SPEC_FILE -out_1D OUTFILE.1D\n"
01842 "\n"
01843 " This program is meant to read in a surface or surface pair,\n"
01844 " and to output and user-requested measures over the surfaces.\n"
01845 " The surfaces must be specified in the SPEC_FILE.\n"
01846 "\n"
01847 " ** Use the 'inspec' command for getting information about the\n"
01848 " surfaces in a spec file.\n"
01849 "\n"
01850 " The output will be a 1D format text file, with one column\n"
01851 " (or possibly 3) per user-specified measure function. Some\n"
01852 " functions require only 1 surface, some require 2.\n"
01853 "\n"
01854 " Current functions (applied with '-func') include:\n"
01855 "\n",
01856 prog, prog );
01857
01858
01859 for ( c = E_SM_INVALID + 1; c < E_SM_FINAL; c++ )
01860 printf( " %-12s : %s\n", g_sm_names[c], g_sm_desc[c]);
01861
01862 printf(
01863 "\n"
01864 "------------------------------------------------------------\n"
01865 "\n"
01866 " examples:\n"
01867 "\n"
01868 " 1. For each node on the surface smoothwm in the spec file,\n"
01869 " fred.spec, output the node number (the default action),\n"
01870 " the xyz coordinates, and the area associated with the\n"
01871 " node (1/3 of the total area of triangles having that node\n"
01872 " as a vertex).\n"
01873 "\n"
01874 " %s \\\n"
01875 " -spec fred1.spec \\\n"
01876 " -sv fred_anat+orig \\\n"
01877 " -surf_A smoothwm \\\n"
01878 " -func coord_A \\\n"
01879 " -func n_area_A \\\n"
01880 " -out_1D fred1_areas.1D \n"
01881 "\n"
01882 " 2. For each node of the surface pair smoothwm and pial,\n"
01883 " display the:\n"
01884 " o node index\n"
01885 " o node's area from the first surface\n"
01886 " o node's area from the second surface\n"
01887 " o node's (approximate) resulting volume\n"
01888 " o thickness at that node (segment distance)\n"
01889 " o coordinates of the first segment node\n"
01890 " o coordinates of the second segment node\n"
01891 "\n"
01892 " Additionally, display total surface areas, minimum and\n"
01893 " maximum thicknesses, and the total volume for the\n"
01894 " cortical ribbon (the sum of node volumes).\n"
01895 "\n"
01896 " %s \\\n"
01897 " -spec fred2.spec \\\n"
01898 " -sv fred_anat+orig \\\n"
01899 " -surf_A smoothwm \\\n"
01900 " -surf_B pial \\\n"
01901 " -func n_area_A \\\n"
01902 " -func n_area_B \\\n"
01903 " -func node_vol \\\n"
01904 " -func thick \\\n"
01905 " -func coord_A \\\n"
01906 " -func coord_B \\\n"
01907 " -info_area \\\n"
01908 " -info_thick \\\n"
01909 " -info_vol \\\n"
01910 " -out_1D fred2_vol.1D \n"
01911 "\n"
01912 " 3. For each node of the surface pair, display the:\n"
01913 " o node index\n"
01914 " o angular diff between the first and second norms\n"
01915 " o angular diff between the segment and first norm\n"
01916 " o angular diff between the segment and second norm\n"
01917 " o the normal vectors for the first surface nodes\n"
01918 " o the normal vectors for the second surface nodes\n"
01919 " o angular diff between the segment and second norm\n"
01920 "\n"
01921 " %s \\\n"
01922 " -spec fred2.spec \\\n"
01923 " -surf_A smoothwm \\\n"
01924 " -surf_B pial \\\n"
01925 " -func ang_norms \\\n"
01926 " -func ang_ns_A \\\n"
01927 " -func ang_ns_B \\\n"
01928 " -func norm_A \\\n"
01929 " -func norm_B \\\n"
01930 " -out_1D fred2_norm_angles.1D \n"
01931 "\n"
01932 " 4. Similar to #3, but output extra debug info, and in\n"
01933 " particular, info regarding node 5000.\n"
01934 "\n"
01935 " %s \\\n"
01936 " -spec fred2.spec \\\n"
01937 " -sv fred_anat+orig \\\n"
01938 " -surf_A smoothwm \\\n"
01939 " -surf_B pial \\\n"
01940 " -func ang_norms \\\n"
01941 " -func ang_ns_A \\\n"
01942 " -func ang_ns_B \\\n"
01943 " -debug 2 \\\n"
01944 " -dnode 5000 \\\n"
01945 " -out_1D fred2_norm_angles.1D \n"
01946 "\n"
01947 " 5. For each node, output the volume, thickness and areas,\n"
01948 " but restrict the nodes to the list contained in column 0\n"
01949 " of file sdata.1D. Furthermore, restrict those nodes to\n"
01950 " the mask inferred by the given '-cmask' option.\n"
01951 "\n"
01952 " %s \\\n"
01953 " -spec fred2.spec \\\n"
01954 " -sv fred_anat+orig \\\n"
01955 " -surf_A smoothwm \\\n"
01956 " -surf_B pial \\\n"
01957 " -func node_vol \\\n"
01958 " -func thick \\\n"
01959 " -func n_area_A \\\n"
01960 " -func n_area_B \\\n"
01961 " -nodes_1D 'sdata.1D[0]' \\\n"
01962 " -cmask '-a sdata.1D[2] -expr step(a-1000)' \\\n"
01963 " -out_1D fred2_masked.1D \n"
01964 "\n"
01965 "------------------------------------------------------------\n"
01966 "\n"
01967 " REQUIRED COMMAND ARGUMENTS:\n"
01968 "\n"
01969 " -spec SPEC_FILE : SUMA spec file\n"
01970 "\n"
01971 " e.g. -spec fred2.spec\n"
01972 "\n"
01973 " The surface specification file contains a list of\n"
01974 " related surfaces. In order for a surface to be\n"
01975 " processed by this program, it must exist in the spec\n"
01976 " file.\n"
01977 "\n"
01978 " -surf_A SURF_NAME : surface name (in spec file)\n"
01979 " -surf_B SURF_NAME : surface name (in spec file)\n"
01980 "\n"
01981 " e.g. -surf_A smoothwm\n"
01982 " e.g. -surf_A lh.smoothwm\n"
01983 " e.g. -surf_B lh.pial\n"
01984 "\n"
01985 " This is used to specify which surface(s) will be used\n"
01986 " by the program. The 'A' and 'B' correspond to other\n"
01987 " program options (e.g. the 'A' in n_area_A).\n"
01988 "\n"
01989 " The '-surf_B' parameter is required only when the user\n"
01990 " wishes to input two surfaces.\n"
01991 "\n"
01992 " Any surface name provided must be unique in the spec\n"
01993 " file, and must match the name of the surface data file\n"
01994 " (e.g. lh.smoothwm.asc).\n"
01995 "\n"
01996 " -out_1D OUT_FILE.1D : 1D output filename\n"
01997 "\n"
01998 " e.g. -out_1D pickle_norm_info.1D\n"
01999 "\n"
02000 " This option is used to specify the name of the output\n"
02001 " file. The output file will be in the 1D ascii format,\n"
02002 " with 2 rows of comments for column headers, and 1 row\n"
02003 " for each node index.\n"
02004 "\n"
02005 " There will be 1 or 3 columns per '-func' option, with\n"
02006 " a default of 1 for \"nodes\".\n"
02007 "\n"
02008 "------------------------------------------------------------\n"
02009 "\n"
02010 " ALPHABETICAL LISTING OF OPTIONS:\n"
02011 "\n"
02012 " -cmask COMMAND : restrict nodes with a mask\n"
02013 "\n"
02014 " e.g. -cmask '-a sdata.1D[2] -expr step(a-1000)'\n"
02015 "\n"
02016 " This option will produce a mask to be applied to the\n"
02017 " list of surface nodes. The total mask size, including\n"
02018 " zero entries, must match the number of nodes. If a\n"
02019 " specific node list is provided via the '-nodes_1D'\n"
02020 " option, then the mask size should match the length of\n"
02021 " the provided node list.\n"
02022 " \n"
02023 " Consider the provided example using the file sdata.1D.\n"
02024 " If a surface has 100000 nodes (and no '-nodes_1D' option\n"
02025 " is used), then there must be 100000 values in column 2\n"
02026 " of the file sdata.1D.\n"
02027 "\n"
02028 " Alternately, if the '-nodes_1D' option is used, giving\n"
02029 " a list of 42 nodes, then the mask length should also be\n"
02030 " 42 (regardless of 0 entries).\n"
02031 "\n"
02032 " See '-nodes_1D' for more information.\n"
02033 "\n"
02034 " -debug LEVEL : display extra run-time info\n"
02035 "\n"
02036 " e.g. -debug 2\n"
02037 " default: -debug 0\n"
02038 "\n"
02039 " Valid debug levels are from 0 to 5.\n"
02040 "\n"
02041 " -dnode NODE : display extra info for node NODE\n"
02042 "\n"
02043 " e.g. -dnode 5000\n"
02044 "\n"
02045 " This option can be used to display extra information\n"
02046 " about node NODE during surface evaluation.\n"
02047 "\n"
02048 " -func FUNCTION : request output for FUNCTION\n"
02049 "\n"
02050 " e.g. -func thick\n"
02051 "\n"
02052 " This option is used to request output for the given\n"
02053 " FUNCTION (measure). Some measures produce one column\n"
02054 " of output (e.g. thick or ang_norms), and some produce\n"
02055 " three (e.g. coord_A). These options, in the order they\n"
02056 " are given, determine the structure of the output file.\n"
02057 "\n"
02058 " Current functions include:\n"
02059 "\n",
02060 prog, prog, prog, prog, prog );
02061
02062
02063 for ( c = E_SM_INVALID + 1; c < E_SM_FINAL; c++ )
02064 printf( " %-12s : %s\n", g_sm_names[c], g_sm_desc[c]);
02065
02066 printf(
02067 "\n"
02068 " Note that the node volumes are approximations. Places\n"
02069 " where either normal points in the 'wrong' direction\n"
02070 " will be incorrect, as will be the parts of the surface\n"
02071 " that 'encompass' this region. Maybe we could refer\n"
02072 " to this as a mushroom effect...\n"
02073 "\n"
02074 " Basically, expect the total volume to be around 10%%\n"
02075 " too large.\n"
02076 "\n"
02077 " ** for more accuracy, try 'SurfPatch -vol' **\n"
02078 "\n"
02079 " -help : show this help menu\n"
02080 "\n"
02081 " -hist : display program revision history\n"
02082 "\n"
02083 " This option is used to provide a history of changes\n"
02084 " to the program, along with version numbers.\n"
02085 "\n"
02086 " NOTE: the following '-info_XXXX' options are used to display\n"
02087 " pieces of 'aggregate' information about the surface(s).\n"
02088 "\n"
02089 " -info_all : display all final info\n"
02090 "\n"
02091 " This is a short-cut to get all '-info_XXXX' options.\n"
02092 "\n"
02093 " -info_area : display info on surface area(s)\n"
02094 "\n"
02095 " Display the total area of each triangulated surface.\n"
02096 "\n"
02097 " -info_norms : display info about the normals\n"
02098 "\n"
02099 " For 1 or 2 surfaces, this will give (if possible) the\n"
02100 " average angular difference between:\n"
02101 "\n"
02102 " o the normals of the surfaces\n"
02103 " o the connecting segment and the first normal\n"
02104 " o the connecting segment and the second normal\n"
02105 "\n"
02106 " -info_thick : display min and max thickness\n"
02107 "\n"
02108 " For 2 surfaces, this is used to display the minimum and\n"
02109 " maximum distances between the surfaces, along each of\n"
02110 " the connecting segments.\n"
02111 "\n"
02112 " -info_vol : display info about the volume\n"
02113 "\n"
02114 " For 2 surfaces, display the total computed volume.\n"
02115 " Note that this node-wise volume computation is an\n"
02116 " approximation, and tends to run ~10 %% high.\n"
02117 "\n"
02118 " ** for more accuracy, try 'SurfPatch -vol' **\n"
02119 "\n"
02120 " -nodes_1D NODELIST.1D : request output for only these nodes\n"
02121 "\n"
02122 " e.g. -nodes_1D node_index_list.1D\n"
02123 " e.g. -nodes_1D sdata.1D'[0]'\n"
02124 "\n"
02125 " The NODELIST file should contain a list of node indices.\n"
02126 " Output from the program would then be restricted to the\n"
02127 " nodes in the list.\n"
02128 " \n"
02129 " For instance, suppose that the file BA_04.1D contains\n"
02130 " a list of surface nodes that are located in Broadman's\n"
02131 " Area 4. To get output from the nodes in that area, use:\n"
02132 " \n"
02133 " -nodes_1D BA_04.1D\n"
02134 " \n"
02135 " For another example, suppose that the file sdata.1D has\n"
02136 " node indices in column 0, and Broadman's Area indices in\n"
02137 " column 3. To restrict output to the nodes in Broadman's\n"
02138 " area 4, use the pair of options:\n"
02139 " \n"
02140 " -nodes_1D 'sdata.1D[0]' \\\n"
02141 " -cmask '-a sdata.1D[3] -expr (1-bool(a-4))' \n"
02142 "\n"
02143 " -sv SURF_VOLUME : specify an associated AFNI volume\n"
02144 "\n"
02145 " e.g. -sv fred_anat+orig\n"
02146 "\n"
02147 " If there is any need to know the orientation of the\n"
02148 " surface, a surface volume dataset may be provided.\n"
02149 "\n"
02150 " -ver : show version information\n"
02151 "\n"
02152 " Show version and compile date.\n"
02153 "\n"
02154 "------------------------------------------------------------\n"
02155 "\n"
02156 " Author: R. Reynolds - %s\n"
02157 "\n",
02158 VERSION );
02159 }
02160 else if ( use_type == ST_USE_HIST )
02161 {
02162 fputs (g_history, stdout);
02163 }
02164 else if ( use_type == ST_USE_VERSION )
02165 {
02166 printf("%s: %s, compile date: %s\n", prog, VERSION, __DATE__);
02167 }
02168 else
02169 fprintf(stderr,"** error: usage - invalid use_type %d\n", use_type);
02170
02171 RETURN(-1);
02172 }
02173
02174
02175
02176
02177
02178
02179 int check_func_name( char * func )
02180 {
02181 int fnum;
02182
02183 ENTRY("check_func_name");
02184
02185 if ( !func )
02186 RETURN(E_SM_INVALID);
02187
02188
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 }
02208
02209
02210
02211
02212
02213
02214 int final_cleanup( opts_t * opts, param_t * p )
02215 {
02216
02217 ENTRY("final_cleanup");
02218
02219
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 }
02257
02258
02259
02260
02261
02262
02263
02264
02265 double tetra_volume( float * p0, float * p1, float * p2, float * p3 )
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 }
02285
02286
02287
02288
02289
02290
02291 double dist_fn( int len, float * p1, float * p2 )
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 }
02312
02313
02314
02315
02316
02317
02318 int disp_surf_t( char * info, surf_t * d )
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 }
02345
02346
02347
02348
02349
02350
02351 int disp_func_t( char * info, func_t * d )
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 }
02371
02372
02373
02374
02375
02376
02377 int disp_opts_t( char * info, opts_t * d )
02378 {
02379 if ( info )
02380 fputs( info, stderr );
02381
02382 if ( ! d )
02383 {
02384 fprintf(stderr,"** disp_opts_t: d == NULL\n");
02385 return -1;
02386 }
02387
02388 fprintf(stderr,
02389 "opts_t struct at %p:\n"
02390 " spec_file = %s\n"
02391 " sv_file = %s\n"
02392 " out_1D_file = %s\n"
02393 " cmask_cmd = %s\n"
02394 " nodes_1D_file = %s\n"
02395 " surf_names[0,1] = %s, %s\n"
02396 " info = %0x\n"
02397 " debug, dnode = %d, %d\n",
02398 d,
02399 CHECK_NULL_STR(d->spec_file),
02400 CHECK_NULL_STR(d->sv_file),
02401 CHECK_NULL_STR(d->out_1D_file),
02402 CHECK_NULL_STR(d->cmask_cmd),
02403 CHECK_NULL_STR(d->nodes_1D_file),
02404 CHECK_NULL_STR(d->surf_names[0]),
02405 CHECK_NULL_STR(d->surf_names[1]),
02406 d->info, d->debug, d->dnode);
02407
02408 return 0;
02409 }
02410
02411
02412
02413
02414
02415
02416 int disp_param_t( char * info, param_t * d )
02417 {
02418 if ( info )
02419 fputs( info, stderr );
02420
02421 if ( ! d )
02422 {
02423 fprintf(stderr,"** disp_param_t: d == NULL\n");
02424 return -1;
02425 }
02426
02427 fprintf(stderr,
02428 "param_t struct at %p:\n"
02429 " F, outfp = %p, %p\n"
02430 " nodes, nnodes = %p, %d\n"
02431 " cmask = %p\n"
02432 " ncmask, ccount = %d, %d\n",
02433 d,
02434 d->F, d->outfp, d->nodes, d->nnodes,
02435 d->cmask, d->ncmask, d->ccount );
02436
02437 return 0;
02438 }
02439
02440
02441
02442
02443
02444
02445 int disp_f3_point( char * info, float * d )
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 }
02462
02463
02464
02465
02466
02467
02468
02469 double lazy_det( int width, double * data )
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
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;
02497
02498 if ( data[c] == 0.0 )
02499 continue;
02500
02501
02502 for ( row = 1; row < width; row++ )
02503 {
02504 rowp = data + row * width;
02505 cur = newm + (row-1) * w2;
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 }
02523
02524
