00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #define _VOL2SURF_C_
00025
00026 char gv2s_history[] =
00027 "----------------------------------------------------------------------\n"
00028 "vol2surf library history:\n"
00029 "\n"
00030 "September 01, 2004 [rickr]\n"
00031 " - initial install into afni\n"
00032 "\n"
00033 "September 02, 2004 [rickr]\n"
00034 " - moved gv2s_map_names here (from SUMA_3dVol2Surf.c)\n"
00035 " - moved v2s_map_type (and test) here (from SUMA_3dVol2Surf.c)\n"
00036 " - define _VOL2SURF_C_, to allow extern defines in vol2surf.h\n"
00037 "\n"
00038 "September 09, 2004 [rickr]\n"
00039 " - in afni_vol2surf(), print v2s options when debug > 1\n"
00040 " - allow (first_node > last_node) if (last == 0), then change to n-1\n"
00041 "\n"
00042 "September 16, 2004 [rickr]\n"
00043 " - added support for -gp_index, computing over a single sub-brick\n"
00044 " - altered subs in dump_surf_3dt(), max_index in set_surf_results(),\n"
00045 " set brick_index in v2s_apply_filter(), mem in alloc_output_mem(),\n"
00046 " and added an index check in validate_v2s_inputs()\n"
00047 " - add gp_index as a parameter of afni_vol2surf()\n"
00048 " - changed keep_norm_dir to norm_dir, allowing default/keep/reverse\n"
00049 "\n"
00050 "September 29, 2004 [rickr]\n"
00051 " - added thd_mask_from_brick() (to make byte mask from sub-brick)\n"
00052 " - added compact_results(), in case nalloc > nused\n"
00053 " - added realloc_ints() and realloc_vals_list() (for compact_results())\n"
00054 " - in afni_vol2surf(), if 1 surf and no norms, set steps to 1\n"
00055 " - in set_surf_results(), pass gp_index to v2s_apply_filter\n"
00056 " - segment oob only if both nodes are\n"
00057 " - move dset bounds to range_3dmm struct\n"
00058 " - in segment_imarr()\n"
00059 " - changed THD_3dmm_to_3dind() to new THD_3dmm_to_3dind_no_wod()\n"
00060 " - verify success of THD_extract_series()\n"
00061 " - keep track of repeated and oob nodes\n"
00062 " - in init_seg_endpoints(), nuke p1, pn; save dicom_to_mm until end\n"
00063 " - changed THD_3dmm_to_3dind() to new THD_3dmm_to_3dind_no_wod()\n"
00064 " - added function v2s_is_good_map_index()\n"
00065 "\n"
00066 "October 08, 2004 [rickr]\n"
00067 " - added disp_v2s_plugin_opts()\n"
00068 " - dealt with default v2s mapping of surface pairs\n"
00069 " - added fill_sopt_default()\n"
00070 " - moved v2s_write_outfile_*() here, with print_header()\n"
00071 " - in afni_vol2surf(), actually write output files\n"
00072 "\n"
00073 "October 25, 2004 [rickr]\n"
00074 " - apply debug and dnode, even for defaults\n"
00075 " - if the user sets dnode, then skip any (debug > 0) tests for it\n"
00076 " - check for out of bounds, even if an endpoint is in (e.g. midpoint)\n"
00077 " - in thd_mask_from_brick(), allow for absolute thresholding\n"
00078 "\n"
00079 "March 22, 2005 [rickr]\n"
00080 " - remove tabs\n"
00081 "---------------------------------------------------------------------\n";
00082
00083 #include "mrilib.h"
00084 #include "vol2surf.h"
00085
00086
00087
00088 typedef struct
00089 {
00090 int nused;
00091 int nalloc;
00092 float * list;
00093 } float_list_t;
00094
00095 typedef struct
00096 {
00097 THD_3dim_dataset * dset;
00098 THD_fvec3 p1;
00099 THD_fvec3 pn;
00100 THD_fvec3 dset_min;
00101 THD_fvec3 dset_max;
00102 int oob_check;
00103 int debug;
00104 } range_3dmm;
00105
00106 typedef struct
00107 {
00108 MRI_IMARR ims;
00109 int repeats;
00110 int masked;
00111 int oob;
00112 int ifirst;
00113 THD_ivec3 i3first;
00114 THD_ivec3 * i3arr;
00115 } range_3dmm_res;
00116
00117
00118
00119 static v2s_results * alloc_output_mem (v2s_opts_t *sopt, v2s_param_t *p);
00120
00121 static int alloc_ints(int ** ptr, int length, char * dstr, int debug);
00122 static int alloc_vals_list(float *** ptr, int length, int width, int debug);
00123 static int check_SUMA_surface( SUMA_surface * s );
00124 static int compact_results(v2s_results * sd, int debug);
00125 static float directed_dist(float * pnew, float * pold, float *dir, float dist);
00126 static float dist_f3mm( THD_fvec3 * p1, THD_fvec3 * p2 );
00127 static int disp_range_3dmm( char * info, range_3dmm * dp );
00128 static int disp_range_3dmm_res( char * info, range_3dmm_res * dp );
00129 static int disp_surf_vals( char * mesg, v2s_results * sd, int node );
00130 static int dump_surf_3dt(v2s_opts_t *sopt, v2s_param_t *p, v2s_results *sd);
00131 static int f3mm_out_of_bounds(THD_fvec3 *cp, THD_fvec3 *min, THD_fvec3 *max);
00132 static int fill_sopt_default(v2s_opts_t * sopt, int nsurf );
00133 static int float_list_alloc(float_list_t *f,int **ilist,int size,int trunc);
00134 static int float_list_comp_mode(float_list_t *f, float *mode, int *nvals,
00135 int *index);
00136 static int float_list_slow_sort(float_list_t * f, int * ilist);
00137 static int init_seg_endpoints(v2s_opts_t * sopt, v2s_param_t * p,
00138 range_3dmm * R, int node );
00139 static int init_range_structs( range_3dmm * r3, range_3dmm_res * res3 );
00140 static double magnitude_f( float * p, int length );
00141 static int print_header(FILE *outfp, char *surf, char *map, v2s_results *sd);
00142 static int realloc_ints( int ** ptr, int length, char * dstr, int debug );
00143 static int realloc_vals_list(float ** ptr, int length, int width, int debug);
00144 static int set_3dmm_bounds(THD_3dim_dataset *dset, THD_fvec3 *min,
00145 THD_fvec3 *max);
00146 static int set_all_surf_vals (v2s_results * sd, int node_ind, int vind,
00147 int i, int j, int k, float fval);
00148 static int set_surf_results(v2s_param_t *p, v2s_opts_t *sopt,v2s_results *sd,
00149 range_3dmm_res * r3res, int node, int findex);
00150 static int segment_imarr(range_3dmm_res *res,range_3dmm *R,v2s_opts_t *sopt,
00151 byte * cmask);
00152 static int v2s_adjust_endpts(v2s_opts_t *sopt, THD_fvec3 *p1, THD_fvec3 *pn);
00153 static float v2s_apply_filter(range_3dmm_res *rr, v2s_opts_t *sopt, int index,
00154 int * findex);
00155 static int v2s_map_needs_sort(int map);
00156 static int validate_v2s_inputs(v2s_opts_t * sopt, v2s_param_t * p);
00157
00158
00159
00160 v2s_plugin_opts gv2s_plug_opts = {0, 0, 0, -1, -1, -1, -1};
00161
00162 char * gv2s_map_names[] = { "none", "mask", "midpoint", "mask2", "ave",
00163 "count", "min", "max", "max_abs", "seg_vals",
00164 "median", "mode" };
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183 v2s_results * afni_vol2surf ( THD_3dim_dataset * gpar, int gp_index,
00184 SUMA_surface * sA, SUMA_surface * sB,
00185 byte * mask, int use_defaults )
00186 {
00187 static v2s_param_t P;
00188 v2s_opts_t * sopt, sopt_def;
00189 v2s_results * res;
00190
00191 ENTRY("afni_vol2surf");
00192
00193 if ( !gpar ) RETURN(NULL);
00194
00195 if ( check_SUMA_surface(sA) ) RETURN(NULL);
00196 if ( sB && check_SUMA_surface(sB) ) RETURN(NULL);
00197
00198 if ( use_defaults )
00199 {
00200 sopt = &sopt_def;
00201 fill_sopt_default(sopt, sB ? 2 : 1);
00202
00203
00204 sopt->debug = gv2s_plug_opts.sopt.debug;
00205 sopt->dnode = gv2s_plug_opts.sopt.dnode;
00206 }
00207 else
00208 sopt = &gv2s_plug_opts.sopt;
00209
00210 sopt->gp_index = gp_index;
00211
00212
00213 memset(&P, 0, sizeof(P));
00214 P.gpar = gpar;
00215 P.cmask = mask;
00216 P.nvox = DSET_NVOX(gpar);
00217 P.over_steps = v2s_vals_over_steps(sopt->map);
00218 P.nsurf = sB ? 2 : 1;
00219 P.surf[0] = *sA;
00220
00221
00222 if ( P.nsurf == 1 && ! sopt->use_norms )
00223 sopt->f_steps = 1;
00224
00225 if ( sB ) P.surf[1] = *sB;
00226
00227 if ( gv2s_plug_opts.sopt.debug > 1 )
00228 disp_v2s_opts_t(" surf options: ", sopt);
00229
00230
00231
00232 res = vol2surf(sopt, &P);
00233
00234
00235 if (res && sopt->outfile_1D )
00236 {
00237 if ( THD_is_file(sopt->outfile_1D) )
00238 fprintf(stderr,"** over-writing 1D output file '%s'\n",
00239 sopt->outfile_1D);
00240 v2s_write_outfile_1D(sopt, res, P.surf[0].label);
00241 }
00242
00243 if (res && sopt->outfile_niml )
00244 {
00245 if ( THD_is_file(sopt->outfile_niml) )
00246 fprintf(stderr,"** over-writing niml output file '%s'\n",
00247 sopt->outfile_niml);
00248 v2s_write_outfile_niml(sopt, res, 0);
00249 }
00250
00251 RETURN(res);
00252 }
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267 v2s_results * vol2surf ( v2s_opts_t * sopt, v2s_param_t * p )
00268 {
00269 v2s_results * sd;
00270 int rv;
00271 ENTRY("vol2surf");
00272
00273 if ( sopt == NULL || p == NULL )
00274 {
00275 fprintf( stderr, "** smd_wo - bad params (%p,%p)\n", sopt, p );
00276 RETURN(NULL);
00277 }
00278
00279 if ( validate_v2s_inputs(sopt, p) )
00280 RETURN(NULL);
00281
00282 if ( sopt->map == E_SMAP_INVALID )
00283 {
00284 fprintf(stderr,"** v2s wo: invalid map %d\n", sopt->map);
00285 RETURN(NULL);
00286 }
00287
00288 sd = alloc_output_mem( sopt, p );
00289 if ( !sd ) RETURN(NULL);
00290
00291 if ( sopt->debug > 1 ) disp_v2s_param_t( "-d post alloc_output_mem : ", p );
00292
00293 rv = dump_surf_3dt( sopt, p, sd );
00294
00295 if ( compact_results(sd, sopt->debug) )
00296 {
00297 free_v2s_results(sd);
00298 RETURN(NULL);
00299 }
00300
00301 if ( sopt->debug > 1 ) disp_v2s_results( "-d post surf creation : ", sd);
00302
00303 RETURN(sd);
00304 }
00305
00306
00307
00308 static int compact_results(v2s_results * sd, int debug)
00309 {
00310 int rv = 0, mem = 0;
00311 ENTRY("compact_results");
00312
00313 if ( sd->nused > sd->nalloc )
00314 {
00315 fprintf(stderr,"** cr: nused (%d) > nalloc (%d) !!\n",
00316 sd->nused, sd->nalloc);
00317 RETURN(-1);
00318 }
00319
00320 if ( sd->nused == sd->nalloc ) RETURN(0);
00321
00322
00323
00324 sd->nalloc = sd->nused;
00325
00326 if ( sd->nodes )
00327 {
00328 mem += sizeof(int);
00329 rv = realloc_ints(&sd->nodes, sd->nalloc, "nodes", debug);
00330 }
00331
00332 if ( ! rv && sd->volind )
00333 {
00334 mem += sizeof(int);
00335 rv = realloc_ints(&sd->volind, sd->nalloc, "volind", debug);
00336 }
00337
00338 if ( ! rv && sd->i )
00339 {
00340 mem += sizeof(int);
00341 rv = realloc_ints(&sd->i, sd->nalloc, "i", debug);
00342 }
00343
00344 if ( ! rv && sd->j )
00345 {
00346 mem += sizeof(int);
00347 rv = realloc_ints(&sd->j, sd->nalloc, "j", debug);
00348 }
00349
00350 if ( ! rv && sd->k )
00351 {
00352 mem += sizeof(int);
00353 rv = realloc_ints(&sd->k, sd->nalloc, "k", debug);
00354 }
00355
00356 if ( ! rv && sd->nvals )
00357 {
00358 mem += sizeof(int);
00359 rv = realloc_ints(&sd->nvals, sd->nalloc, "nvals", debug);
00360 }
00361
00362 if ( ! rv )
00363 {
00364 mem += (sizeof(float) * sd->max_vals);
00365 rv = realloc_vals_list(sd->vals, sd->nalloc, sd->max_vals, debug);
00366 }
00367
00368 if ( rv ) RETURN(rv);
00369
00370 mem *= sd->nalloc;
00371
00372 if ( debug > 1 )
00373 fprintf(stderr,"+d compact results: reallocated %d bytes down to %d\n",
00374 sd->memory, mem);
00375
00376 sd->memory = mem;
00377
00378 RETURN(rv);
00379 }
00380
00381
00382
00383
00384
00385
00386
00387 static int dump_surf_3dt( v2s_opts_t * sopt, v2s_param_t * p, v2s_results * sd )
00388 {
00389 range_3dmm_res r3mm_res;
00390 range_3dmm r3mm;
00391 float dist, min_dist, max_dist;
00392 int sub, subs, nindex, findex = 0;
00393 int oobc, oomc, max_index;
00394 int oob1, oob2;
00395
00396 ENTRY("dump_surf_3dt");
00397
00398 if ( ! sopt || ! p || ! sd )
00399 {
00400 fprintf(stderr, "** ds3 : bad params (%p,%p,%p)\n", sopt, p, sd );
00401 RETURN(-1);
00402 }
00403
00404
00405 subs = sopt->gp_index >= 0 ? 1 : DSET_NVALS(p->gpar);
00406 init_range_structs( &r3mm, &r3mm_res );
00407 r3mm.dset = p->gpar;
00408 set_3dmm_bounds( p->gpar, &r3mm.dset_min, &r3mm.dset_max );
00409
00410 if ( sopt->debug > 1 )
00411 fprintf(stderr, "-d dset bounding box: (%f, %f, %f)\n"
00412 " (%f, %f, %f)\n",
00413 r3mm.dset_min.xyz[0],r3mm.dset_min.xyz[1],r3mm.dset_min.xyz[2],
00414 r3mm.dset_max.xyz[0],r3mm.dset_max.xyz[1],r3mm.dset_max.xyz[2]);
00415
00416 min_dist = 9999.9;
00417 max_dist = -1.0;
00418 oobc = 0;
00419 oomc = 0;
00420
00421
00422
00423 for ( nindex = sopt->first_node; nindex <= sopt->last_node; nindex++ )
00424 {
00425
00426 max_index = p->over_steps ? sopt->f_steps : subs;
00427
00428 init_seg_endpoints(sopt, p, &r3mm, nindex);
00429 v2s_adjust_endpts( sopt, &r3mm.p1, &r3mm.pn );
00430
00431 if ( r3mm.debug )
00432 r3mm.debug = 0;
00433
00434 if ( nindex == sopt->dnode )
00435 r3mm.debug = sopt->debug > 0 ? sopt->debug : 1;
00436
00437
00438 oob1 = f3mm_out_of_bounds( &r3mm.p1, &r3mm.dset_min, &r3mm.dset_max );
00439 oob2 = f3mm_out_of_bounds( &r3mm.pn, &r3mm.dset_min, &r3mm.dset_max );
00440 if ( oob1 && oob2 )
00441 {
00442 oobc++;
00443 if ( sopt->oob.show )
00444 if ( set_all_surf_vals( sd, nindex, sopt->oob.index,
00445 sopt->oob.index, sopt->oob.index,
00446 sopt->oob.index, sopt->oob.value) )
00447 RETURN(1);
00448 if ( nindex == sopt->dnode )
00449 {
00450 disp_surf_vals("-d debug node, out-of-bounds : ", sd, -1);
00451 fprintf(stderr,"-d dnode coords: (%f, %f, %f)\n",
00452 r3mm.p1.xyz[0], r3mm.p1.xyz[1], r3mm.p1.xyz[2]);
00453 }
00454 continue;
00455 }
00456 else
00457 r3mm.oob_check = oob1 || oob2;
00458
00459 dist = dist_f3mm( &r3mm.p1, &r3mm.pn );
00460 if ( dist < min_dist ) min_dist = dist;
00461 if ( dist > max_dist ) max_dist = dist;
00462
00463 if ( segment_imarr( &r3mm_res, &r3mm, sopt, p->cmask ) != 0 )
00464 continue;
00465
00466 if ( r3mm_res.ims.num == 0 )
00467 {
00468
00469 if ( r3mm_res.oob == sopt->f_steps )
00470 {
00471 oobc++;
00472 if ( sopt->oob.show )
00473 if ( set_all_surf_vals( sd, nindex, sopt->oob.index,
00474 sopt->oob.index, sopt->oob.index,
00475 sopt->oob.index, sopt->oob.value) )
00476 RETURN(1);
00477 if ( nindex == sopt->dnode )
00478 disp_surf_vals("-d debug node, out-of-bounds : ", sd, -1);
00479 }
00480 else
00481 {
00482 oomc++;
00483 if ( sopt->oom.show )
00484 if ( set_all_surf_vals( sd, nindex, r3mm_res.ifirst,
00485 r3mm_res.i3first.ijk[0], r3mm_res.i3first.ijk[1],
00486 r3mm_res.i3first.ijk[2], sopt->oom.value ) )
00487 RETURN(1);
00488 if ( nindex == sopt->dnode )
00489 disp_surf_vals("-d debug node, out-of-mask : ", sd, -1);
00490 }
00491
00492 if ( nindex == sopt->dnode )
00493 fprintf(stderr,"-d dnode coords: (%f, %f, %f)\n",
00494 r3mm.p1.xyz[0], r3mm.p1.xyz[1], r3mm.p1.xyz[2]);
00495
00496 continue;
00497 }
00498
00499
00500 (void)v2s_apply_filter(&r3mm_res, sopt, 0, &findex);
00501
00502 if ( set_surf_results(p, sopt, sd, &r3mm_res, nindex, findex) )
00503 RETURN(-1);
00504
00505
00506 for ( sub = 0; sub < r3mm_res.ims.num; sub++ )
00507 {
00508 mri_free(r3mm_res.ims.imarr[sub]);
00509 r3mm_res.ims.imarr[sub] = NULL;
00510 }
00511 r3mm_res.ims.num = 0;
00512 }
00513
00514 if ( sopt->debug > 0 )
00515 {
00516 fprintf( stderr, "-d node pair dist (min,max) = (%f,%f)\n",
00517 min_dist, max_dist );
00518 fprintf( stderr, "-d out-of-bounds, o-o-mask counts : %d, %d (of %d)\n",
00519 oobc, oomc, sopt->last_node - sopt->first_node + 1);
00520 }
00521
00522
00523 if ( r3mm_res.ims.nall > 0 )
00524 {
00525 free(r3mm_res.ims.imarr);
00526 free(r3mm_res.i3arr);
00527 r3mm_res.ims.nall = 0;
00528 }
00529
00530 RETURN(0);
00531 }
00532
00533
00534
00535
00536
00537
00538 static int set_surf_results(v2s_param_t *p, v2s_opts_t * sopt, v2s_results * sd,
00539 range_3dmm_res * r3res, int node, int findex)
00540 {
00541 int i, j, k, volind;
00542 int c, max_index;
00543 ENTRY("set_surf_results");
00544
00545 if ( sd->nused >= sd->nalloc )
00546 {
00547 fprintf(stderr,"** ssr: nused (%d) >= nalloc (%d)!\n",
00548 sd->nused, sd->nalloc);
00549 RETURN(1);
00550 }
00551
00552 i = r3res->i3arr[findex].ijk[0];
00553 j = r3res->i3arr[findex].ijk[1];
00554 k = r3res->i3arr[findex].ijk[2];
00555
00556
00557 volind = i + DSET_NX(p->gpar) * (j + DSET_NY(p->gpar) * k );
00558
00559
00560 if (sd->nodes ) sd->nodes[sd->nused] = node;
00561 if (sd->volind) sd->volind[sd->nused] = volind;
00562 if (sd->i ) sd->i[sd->nused] = i;
00563 if (sd->j ) sd->j[sd->nused] = j;
00564 if (sd->k ) sd->k[sd->nused] = k;
00565 if (sd->nvals ) sd->nvals[sd->nused] = r3res->ims.num;
00566
00567
00568 max_index = p->over_steps ? r3res->ims.num : DSET_NVALS(p->gpar);
00569 if ( max_index > sd->max_vals ) max_index = sd->max_vals;
00570
00571 if ( sopt->gp_index >= 0 )
00572 sd->vals[0][sd->nused] =
00573 v2s_apply_filter(r3res, sopt, sopt->gp_index, NULL);
00574 else
00575 for ( c = 0; c < max_index; c++ )
00576 sd->vals[c][sd->nused] = v2s_apply_filter(r3res, sopt, c, NULL);
00577
00578
00579 if ( max_index < sd->max_vals )
00580 for ( c = max_index; c < sd->max_vals; c++ )
00581 sd->vals[c][sd->nused] = 0.0;
00582
00583
00584 if ( node == sopt->dnode )
00585 {
00586 fprintf(stderr,
00587 "--------------------------------------------------\n");
00588 if ( ! p->over_steps && sopt->gp_index >= 0 )
00589 fprintf(stderr,"+d dnode %d gets %f from gp_index %d\n",
00590 node, sd->vals[0][sd->nused], sopt->gp_index);
00591 if ( sopt->debug > 1 )
00592 fprintf(stderr, "-d debug: node %d, findex %d, vol_index %d\n",
00593 node, findex, volind );
00594 if ( sopt->use_norms )
00595 {
00596 float * fp = p->surf[0].norm[node].xyz;
00597 fprintf(stderr,"-d normal %f, %f, %f\n", fp[0], fp[1], fp[2]);
00598 }
00599 disp_mri_imarr( "+d raw data: ", &r3res->ims );
00600 }
00601
00602 sd->nused++;
00603
00604 RETURN(0);
00605 }
00606
00607
00608
00609
00610
00611
00612
00613 static int set_all_surf_vals( v2s_results * sd, int node_ind, int vind,
00614 int i, int j, int k, float fval )
00615 {
00616 int c, nused;
00617
00618 ENTRY("set_all_surf_vals");
00619
00620 nused = sd->nused;
00621
00622 if ( nused >= sd->nalloc )
00623 {
00624 fprintf(stderr,"** sasv: nused=%d >= nalloc=%d!\n", nused, sd->nalloc);
00625 RETURN(1);
00626 }
00627
00628 if ( sd->nodes ) sd->nodes[nused] = node_ind;
00629 if ( sd->volind ) sd->volind[nused] = vind;
00630 if ( sd->i ) sd->i[nused] = i;
00631 if ( sd->j ) sd->j[nused] = j;
00632 if ( sd->k ) sd->k[nused] = k;
00633 if ( sd->nvals ) sd->nvals[nused] = sd->max_vals;
00634
00635 for ( c = 0; c < sd->max_vals; c++ )
00636 sd->vals[c][nused] = fval;
00637
00638 sd->nused++;
00639
00640 RETURN(0);
00641 }
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654 static int segment_imarr( range_3dmm_res *res, range_3dmm *R, v2s_opts_t *sopt,
00655 byte * cmask )
00656 {
00657 THD_fvec3 f3mm;
00658 THD_ivec3 i3ind;
00659 float rat1, ratn;
00660 int nx, ny;
00661 int step, vindex, prev_ind;
00662
00663 ENTRY("segment_imarr");
00664
00665
00666 if ( !R || !sopt || !res || !R->dset )
00667 {
00668 fprintf(stderr, "** seg_imarr: invalid params (%p,%p,%p)\n",R,sopt,res);
00669 if ( R ) disp_range_3dmm("segment_imarr: bad inputs:", R );
00670 RETURN(-1);
00671 }
00672
00673 if ( R->debug > 1 )
00674 disp_range_3dmm("-d segment_imarr: ", R );
00675
00676
00677 if ( sopt->f_steps < 1 )
00678 {
00679 res->ims.num = 0;
00680 RETURN(0);
00681 }
00682
00683 nx = DSET_NX(R->dset);
00684 ny = DSET_NY(R->dset);
00685
00686
00687 if ( res->ims.nall < sopt->f_steps )
00688 {
00689 if ( R->debug > 1 )
00690 fprintf(stderr,"+d realloc of imarr (from %d to %d pointers)\n",
00691 res->ims.nall,sopt->f_steps);
00692
00693 res->ims.nall = sopt->f_steps;
00694 res->ims.imarr = realloc(res->ims.imarr,
00695 sopt->f_steps*sizeof(MRI_IMAGE *));
00696 res->i3arr = realloc(res->i3arr, sopt->f_steps*sizeof(THD_ivec3));
00697 if ( !res->ims.imarr || !res->i3arr )
00698 {
00699 fprintf(stderr,"** seg_imarr: no memory for %d MRI_IMAGE ptrs\n",
00700 sopt->f_steps);
00701 res->ims.nall = res->ims.num = 0;
00702
00703 if ( res->ims.imarr ) free(res->ims.imarr);
00704 if ( res->i3arr ) free(res->i3arr);
00705 RETURN(-1);
00706 }
00707 }
00708
00709
00710 res->ims.num = 0;
00711 res->repeats = 0;
00712 res->masked = 0;
00713 res->oob = 0;
00714 res->i3first = THD_3dmm_to_3dind_no_wod( R->dset, R->p1 );
00715 res->ifirst = res->i3first.ijk[0] +
00716 nx * (res->i3first.ijk[1] + ny * res->i3first.ijk[2] );
00717
00718 prev_ind = -1;
00719
00720 for ( step = 0; step < sopt->f_steps; step++ )
00721 {
00722
00723 if ( sopt->f_steps < 2 )
00724 ratn = 0.0;
00725 else
00726 ratn = (float)step / (sopt->f_steps - 1);
00727 rat1 = 1.0 - ratn;
00728
00729 f3mm.xyz[0] = rat1 * R->p1.xyz[0] + ratn * R->pn.xyz[0];
00730 f3mm.xyz[1] = rat1 * R->p1.xyz[1] + ratn * R->pn.xyz[1];
00731 f3mm.xyz[2] = rat1 * R->p1.xyz[2] + ratn * R->pn.xyz[2];
00732
00733
00734 if ( R->oob_check &&
00735 f3mm_out_of_bounds( &f3mm, &R->dset_min, &R->dset_max ) )
00736 {
00737 res->oob++;
00738 continue;
00739 }
00740
00741 i3ind = THD_3dmm_to_3dind_no_wod( R->dset, f3mm );
00742 vindex = i3ind.ijk[0] + nx * (i3ind.ijk[1] + ny * i3ind.ijk[2] );
00743
00744
00745 if ( cmask && !cmask[vindex] )
00746 {
00747 res->masked++;
00748 continue;
00749 }
00750
00751
00752 if ( sopt->f_index == V2S_INDEX_VOXEL && vindex == prev_ind )
00753 {
00754 res->repeats++;
00755 continue;
00756 }
00757
00758
00759
00760 prev_ind = vindex;
00761
00762
00763
00764 res->i3arr [res->ims.num] = i3ind;
00765 res->ims.imarr[res->ims.num] = THD_extract_series( vindex, R->dset, 0 );
00766 res->ims.num++;
00767
00768 if ( ! res->ims.imarr[res->ims.num-1] )
00769 {
00770 fprintf(stderr,"** failed THD_extract_series, vox %d\n", vindex);
00771 RETURN(-1);
00772 }
00773
00774 if ( R->debug > 2 )
00775 fprintf(stderr, "-d seg step %2d, vindex %d, coords %f %f %f\n",
00776 step,vindex,f3mm.xyz[0],f3mm.xyz[1],f3mm.xyz[2]);
00777 }
00778
00779 if ( R->debug > 1 )
00780 disp_range_3dmm_res( "+d i3mm_seg_imarr results: ", res );
00781
00782 RETURN(0);
00783 }
00784
00785
00786
00787
00788
00789
00790
00791
00792 static int f3mm_out_of_bounds( THD_fvec3 * cp, THD_fvec3 * min, THD_fvec3 * max)
00793 {
00794 int c;
00795
00796 if ( !cp || !min || !max )
00797 return(-1);
00798
00799 for ( c = 0; c < 3; c++ )
00800 {
00801 if ( ( cp->xyz[c] < min->xyz[c] ) ||
00802 ( cp->xyz[c] > max->xyz[c] ) )
00803 return(-1);
00804 }
00805
00806 return(0);
00807 }
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817 static int v2s_adjust_endpts( v2s_opts_t * sopt, THD_fvec3 * p1, THD_fvec3 * pn)
00818 {
00819 THD_fvec3 f3_diff;
00820 float dist, factor;
00821
00822 ENTRY("v2s_adjust_endpts");
00823
00824 if ( !sopt || !p1 || !pn )
00825 {
00826 fprintf(stderr,"** v2s_ae: invalid params (%p,%p,%p)\n", sopt, p1, pn);
00827 RETURN(-1);
00828 }
00829
00830
00831 f3_diff.xyz[0] = pn->xyz[0] - p1->xyz[0];
00832 f3_diff.xyz[1] = pn->xyz[1] - p1->xyz[1];
00833 f3_diff.xyz[2] = pn->xyz[2] - p1->xyz[2];
00834
00835 dist = dist_f3mm( p1, pn );
00836
00837 if ( (sopt->f_p1_fr != 0.0) || (sopt->f_p1_mm != 0.0) )
00838 {
00839 if ( sopt->f_p1_fr != 0.0 )
00840 factor = sopt->f_p1_fr;
00841 else
00842 factor = (dist == 0.0) ? 0.0 : sopt->f_p1_mm / dist;
00843
00844 p1->xyz[0] += factor * f3_diff.xyz[0];
00845 p1->xyz[1] += factor * f3_diff.xyz[1];
00846 p1->xyz[2] += factor * f3_diff.xyz[2];
00847 }
00848
00849 if ( (sopt->f_pn_fr != 0.0) || (sopt->f_pn_mm != 0.0) )
00850 {
00851 if ( sopt->f_pn_fr != 0.0 )
00852 factor = sopt->f_pn_fr;
00853 else
00854 factor = (dist == 0.0) ? 0.0 : sopt->f_pn_mm / dist;
00855
00856 pn->xyz[0] += factor * f3_diff.xyz[0];
00857 pn->xyz[1] += factor * f3_diff.xyz[1];
00858 pn->xyz[2] += factor * f3_diff.xyz[2];
00859 }
00860
00861 switch ( sopt->map )
00862 {
00863 default:
00864 fprintf(stderr,"** v2s_ae: mapping %d not ready\n", sopt->map );
00865 RETURN(-1);
00866
00867 case E_SMAP_AVE:
00868 case E_SMAP_MAX:
00869 case E_SMAP_MAX_ABS:
00870 case E_SMAP_MIN:
00871 case E_SMAP_MASK:
00872 case E_SMAP_SEG_VALS:
00873 case E_SMAP_MEDIAN:
00874 case E_SMAP_MODE:
00875 break;
00876
00877 case E_SMAP_MIDPT:
00878
00879
00880 p1->xyz[0] = (p1->xyz[0] + pn->xyz[0]) / 2.0;
00881 p1->xyz[1] = (p1->xyz[1] + pn->xyz[1]) / 2.0;
00882 p1->xyz[2] = (p1->xyz[2] + pn->xyz[2]) / 2.0;
00883
00884 break;
00885 }
00886
00887 RETURN(0);
00888 }
00889
00890
00891
00892
00893
00894
00895
00896
00897 static float v2s_apply_filter( range_3dmm_res * rr, v2s_opts_t * sopt,
00898 int index, int * findex )
00899 {
00900 static float_list_t flist = { 0, 0, NULL };
00901 static int * ind_list = NULL;
00902 double tmp, comp = 0.0;
00903 float fval;
00904 int count, source;
00905 int brick_index = 0;
00906
00907 ENTRY("v2s_apply_filter");
00908
00909 if ( !rr || !sopt || index < 0 )
00910 {
00911 fprintf(stderr,"** v2s_cm2: invalid params (%p,%p,%d)\n",
00912 rr, sopt, index);
00913 RETURN(0.0);
00914 }
00915
00916 if ( rr->ims.num <= 0 )
00917 RETURN(0.0);
00918
00919
00920 if ( v2s_map_needs_sort( sopt->map ) )
00921 {
00922 if ( float_list_alloc( &flist, &ind_list, rr->ims.num, 0 ) != 0 )
00923 RETURN(0.0);
00924 for ( count = 0; count < rr->ims.num; count++ )
00925 {
00926 flist.list[count] = MRI_FLOAT_PTR(rr->ims.imarr[count])[index];
00927 ind_list [count] = count;
00928 }
00929 flist.nused = rr->ims.num;
00930 float_list_slow_sort( &flist, ind_list );
00931 }
00932
00933 switch ( sopt->map )
00934 {
00935 default:
00936 if ( findex ) *findex = 0;
00937 RETURN(0.0);
00938
00939 case E_SMAP_AVE:
00940 if ( findex ) *findex = 0;
00941 for ( count = 0; count < rr->ims.num; count++ )
00942 comp += MRI_FLOAT_PTR(rr->ims.imarr[count])[index];
00943
00944 comp = comp / rr->ims.num;
00945 break;
00946
00947 case E_SMAP_MASK:
00948 case E_SMAP_MIDPT:
00949 if ( findex ) *findex = 0;
00950
00951 comp = MRI_FLOAT_PTR(rr->ims.imarr[0])[index];
00952 break;
00953
00954 case E_SMAP_MAX:
00955 comp = MRI_FLOAT_PTR(rr->ims.imarr[0])[index];
00956 if ( findex ) *findex = 0;
00957
00958 for ( count = 1; count < rr->ims.num; count++ )
00959 {
00960 tmp = MRI_FLOAT_PTR(rr->ims.imarr[count])[index];
00961 if ( tmp > comp )
00962 {
00963 if ( findex ) *findex = count;
00964 comp = tmp;
00965 }
00966 }
00967 break;
00968
00969 case E_SMAP_MAX_ABS:
00970 comp = MRI_FLOAT_PTR(rr->ims.imarr[0])[index];
00971 if ( findex ) *findex = 0;
00972
00973 for ( count = 1; count < rr->ims.num; count++ )
00974 {
00975 tmp = MRI_FLOAT_PTR(rr->ims.imarr[count])[index];
00976 if ( fabs(tmp) > fabs(comp) )
00977 {
00978 if ( findex ) *findex = count;
00979 comp = tmp;
00980 }
00981 }
00982 break;
00983
00984 case E_SMAP_MIN:
00985 comp = MRI_FLOAT_PTR(rr->ims.imarr[0])[index];
00986 if ( findex ) *findex = 0;
00987
00988 for ( count = 1; count < rr->ims.num; count++ )
00989 {
00990 tmp = MRI_FLOAT_PTR(rr->ims.imarr[count])[index];
00991 if ( tmp < comp )
00992 {
00993 if ( findex ) *findex = count;
00994 comp = tmp;
00995 }
00996 }
00997 break;
00998
00999 case E_SMAP_SEG_VALS:
01000
01001 if ( sopt->gp_index >= 0 ) brick_index = sopt->gp_index;
01002 if ( findex ) *findex = 0;
01003 comp = MRI_FLOAT_PTR(rr->ims.imarr[index])[brick_index];
01004 break;
01005
01006 case E_SMAP_MEDIAN:
01007 count = flist.nused >> 1;
01008 if ( (flist.nused & 1) || (count == 0) )
01009 {
01010 comp = flist.list[count];
01011 if ( findex ) *findex = ind_list[count];
01012 }
01013 else
01014 {
01015 comp = (flist.list[count-1] + flist.list[count]) / 2;
01016 if ( findex ) *findex = ind_list[count-1];
01017 }
01018 break;
01019
01020 case E_SMAP_MODE:
01021 float_list_comp_mode(&flist, &fval, &count, &source);
01022 comp = fval;
01023 if ( findex ) *findex = ind_list[source];
01024 break;
01025 }
01026
01027 RETURN((float)comp);
01028 }
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038 static int v2s_map_needs_sort( int map )
01039 {
01040 if ( (map == E_SMAP_MEDIAN) || (map == E_SMAP_MODE) )
01041 return 1;
01042
01043 return 0;
01044 }
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054 static int float_list_comp_mode( float_list_t *f, float *mode, int *nvals,
01055 int *index )
01056 {
01057 float fcur;
01058 int ncur, c;
01059
01060 ENTRY("float_list_comp_mode");
01061
01062
01063 *nvals = ncur = 1;
01064 *mode = fcur = f->list[0];
01065 *index = 0;
01066
01067 for ( c = 1; c < f->nused; c++ )
01068 {
01069 if ( f->list[c] == fcur )
01070 ncur ++;
01071 else
01072 {
01073 if ( ncur > *nvals )
01074 {
01075 *mode = fcur;
01076 *nvals = ncur;
01077 *index = c;
01078 }
01079
01080 fcur = f->list[c];
01081 ncur = 1;
01082 }
01083 }
01084
01085 if ( ncur > *nvals )
01086 {
01087 *mode = fcur;
01088 *nvals = ncur;
01089 *index = c;
01090 }
01091
01092 RETURN(0);
01093 }
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105 static int float_list_slow_sort( float_list_t * f, int * ilist )
01106 {
01107 float * list, save;
01108 int c0, c1, sindex;
01109
01110 ENTRY("float_list_slow_sort");
01111
01112 list = f->list;
01113
01114 for ( c0 = 0; c0 < f->nused-1; c0++ )
01115 {
01116 sindex = c0;
01117 save = list[c0];
01118
01119
01120 for ( c1 = c0+1; c1 < f->nused; c1++ )
01121 if ( list[c1] < save )
01122 {
01123 sindex = c1;
01124 save = list[sindex];
01125 }
01126
01127
01128 if ( sindex > c0 )
01129 {
01130 list[sindex] = list[c0];
01131 list[c0] = save;
01132
01133 if ( ilist )
01134 {
01135 c1 = ilist[sindex];
01136 ilist[sindex] = ilist[c0];
01137 ilist[c0] = c1;
01138 }
01139 }
01140 }
01141
01142 RETURN(0);
01143 }
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156 static int float_list_alloc(float_list_t *f, int **ilist, int size, int trunc)
01157 {
01158 ENTRY("float_list_alloc");
01159
01160 if ( (f->nalloc < size) ||
01161 (trunc && (f->nalloc > size)) )
01162 {
01163 f->list = (float *)realloc(f->list, size * sizeof(float));
01164 if ( ! f->list )
01165 {
01166 fprintf(stderr,"** float_list_alloc: failed for %d floats\n", size);
01167 RETURN(-2);
01168 }
01169 f->nalloc = size;
01170
01171 if ( ilist )
01172 {
01173 *ilist = (int *)realloc(*ilist, size * sizeof(int));
01174 if ( ! *ilist )
01175 {
01176 fprintf(stderr,"** float_list_alloc: failed for %d ints\n",
01177 size);
01178 RETURN(-2);
01179 }
01180 }
01181
01182 if ( trunc && (f->nused > f->nalloc) )
01183 f->nused = f->nalloc;
01184 }
01185
01186 RETURN(0);
01187 }
01188
01189
01190
01191
01192
01193
01194 static float directed_dist(float * pnew, float * pold, float * dir, float dist)
01195 {
01196 double mag, ratio;
01197 int c;
01198
01199 ENTRY("directed_dist");
01200
01201 mag = magnitude_f(dir, 3);
01202
01203 if ( mag < V2S_EPSILON )
01204 ratio = 0.0;
01205 else
01206 ratio = dist/mag;
01207
01208 for ( c = 0; c < 3; c++ )
01209 pnew[c] = pold[c] + dir[c]*ratio;
01210
01211 RETURN(dist);
01212 }
01213
01214
01215
01216
01217
01218
01219 static int init_seg_endpoints ( v2s_opts_t * sopt, v2s_param_t * p,
01220 range_3dmm * R, int node )
01221 {
01222 SUMA_surface * sp;
01223 ENTRY("init_seg_endpoints");
01224
01225
01226 sp = p->surf;
01227 R->p1.xyz[0] = sp->ixyz[node].x;
01228 R->p1.xyz[1] = sp->ixyz[node].y;
01229 R->p1.xyz[2] = sp->ixyz[node].z;
01230
01231
01232 if ( sopt->use_norms )
01233 directed_dist(R->pn.xyz, R->p1.xyz, sp->norm[node].xyz, sopt->norm_len);
01234 else if ( p->nsurf > 1 )
01235 {
01236
01237 sp = p->surf + 1;
01238 R->pn.xyz[0] = sp->ixyz[node].x;
01239 R->pn.xyz[1] = sp->ixyz[node].y;
01240 R->pn.xyz[2] = sp->ixyz[node].z;
01241 }
01242 else
01243 R->pn = R->p1;
01244
01245 R->p1 = THD_dicomm_to_3dmm(R->dset, R->p1);
01246 R->pn = THD_dicomm_to_3dmm(R->dset, R->pn);
01247
01248 RETURN(0);
01249 }
01250
01251
01252
01253
01254
01255
01256 static int init_range_structs( range_3dmm * r3, range_3dmm_res * res3 )
01257 {
01258 ENTRY("init_range_structs");
01259
01260 r3->dset = NULL;
01261 r3->debug = 0;
01262
01263 res3->ims.num = 0;
01264 res3->ims.nall = 0;
01265 res3->ims.imarr = NULL;
01266 res3->i3arr = NULL;
01267
01268 RETURN(0);
01269 }
01270
01271
01272
01273
01274
01275
01276
01277
01278 int set_3dmm_bounds ( THD_3dim_dataset *dset, THD_fvec3 *min, THD_fvec3 *max)
01279 {
01280 float tmp;
01281 int c;
01282
01283 ENTRY("set_3dmm_bounds");
01284
01285 if ( !dset || !min || !max )
01286 {
01287 fprintf(stderr, "** invalid params to set_3dmm_bounds: (%p,%p,%p)\n",
01288 dset, min, max );
01289 RETURN(-1);
01290 }
01291
01292
01293 min->xyz[0] = DSET_XORG(dset) - 0.5 * DSET_DX(dset);
01294 max->xyz[0] = min->xyz[0] + DSET_NX(dset) * DSET_DX(dset);
01295
01296 min->xyz[1] = DSET_YORG(dset) - 0.5 * DSET_DY(dset);
01297 max->xyz[1] = min->xyz[1] + DSET_NY(dset) * DSET_DY(dset);
01298
01299 min->xyz[2] = DSET_ZORG(dset) - 0.5 * DSET_DZ(dset);
01300 max->xyz[2] = min->xyz[2] + DSET_NZ(dset) * DSET_DZ(dset);
01301
01302 for ( c = 0; c < 3; c++ )
01303 if ( min->xyz[c] > max->xyz[c] )
01304 {
01305 tmp = min->xyz[c];
01306 min->xyz[c] = max->xyz[c];
01307 max->xyz[c] = tmp;
01308 }
01309
01310 RETURN(0);
01311 }
01312
01313
01314
01315
01316
01317
01318 static v2s_results * alloc_output_mem(v2s_opts_t *sopt, v2s_param_t *p)
01319 {
01320 v2s_results * sd;
01321 int rv = 0, mem, nnodes;
01322
01323 ENTRY("allocate_output_mem");
01324
01325
01326 if ( sopt->first_node > sopt->last_node && sopt->last_node > 0 )
01327 {
01328 fprintf(stderr,"** error : first_node (%d) > last_node (%d)\n",
01329 sopt->first_node, sopt->last_node);
01330 RETURN(NULL);
01331 }
01332
01333 nnodes = p->surf[0].num_ixyz;
01334
01335 sd = calloc(1, sizeof(v2s_results));
01336 if ( ! sd )
01337 {
01338 fprintf(stderr,"** aom: failed to allocate v2s_results struct\n");
01339 RETURN(NULL);
01340 }
01341
01342
01343 sd->nodes = sd->volind = sd->i = sd->j = sd->k = sd->nvals = NULL;
01344 sd->vals = NULL;
01345
01346
01347 if ( sopt->first_node < 0 ) sopt->first_node = 0;
01348 if ( sopt->first_node >= nnodes ) sopt->first_node = nnodes-1;
01349 if ( sopt->last_node <= 0 ) sopt->last_node = nnodes-1;
01350 if ( sopt->last_node >= nnodes ) sopt->last_node = nnodes-1;
01351
01352 sd->nalloc = sopt->last_node - sopt->first_node + 1;
01353 sd->nused = 0;
01354
01355
01356 if ( p->over_steps ) sd->max_vals = sopt->f_steps;
01357 else if ( sopt->gp_index >= 0 ) sd->max_vals = 1;
01358 else sd->max_vals = DSET_NVALS(p->gpar);
01359
01360 if ( sopt->skip_cols & V2S_SKIP_VALS ) sd->max_vals = 1;
01361
01362 if ( p->over_steps && (sopt->skip_cols & V2S_SKIP_VALS) )
01363 {
01364 fprintf(stderr,"** if SKIP_VALS, cannot compute results over steps\n");
01365 free(sd);
01366 RETURN(NULL);
01367 }
01368
01369 sd->memory = 0;
01370
01371
01372 mem = 0;
01373 if ( !(sopt->skip_cols & V2S_SKIP_NODES ) ) mem += sizeof(int);
01374 if ( !(sopt->skip_cols & V2S_SKIP_VOLIND) ) mem += sizeof(int);
01375 if ( !(sopt->skip_cols & V2S_SKIP_I ) ) mem += sizeof(int);
01376 if ( !(sopt->skip_cols & V2S_SKIP_J ) ) mem += sizeof(int);
01377 if ( !(sopt->skip_cols & V2S_SKIP_K ) ) mem += sizeof(int);
01378 if ( !(sopt->skip_cols & V2S_SKIP_NVALS ) ) mem += sizeof(int);
01379
01380
01381 if ( !(sopt->skip_cols & V2S_SKIP_VALS ) ) mem += sizeof(float);
01382 else mem += sd->max_vals * sizeof(float);
01383
01384 mem *= sd->nalloc;
01385 sd->memory = mem;
01386
01387 if ( sopt->debug > 1 )
01388 fprintf(stderr,"+d alloc result mem: %d bytes, max_vals = %d\n",
01389 mem, sd->max_vals);
01390
01391
01392
01393 if ( ! (sopt->skip_cols & V2S_SKIP_NODES) )
01394 rv = alloc_ints(&sd->nodes, sd->nalloc, "nodes", sopt->debug);
01395
01396 if ( ! rv && ! (sopt->skip_cols & V2S_SKIP_VOLIND) )
01397 rv = alloc_ints(&sd->volind, sd->nalloc, "volind", sopt->debug);
01398
01399 if ( ! rv && ! (sopt->skip_cols & V2S_SKIP_I) )
01400 rv = alloc_ints(&sd->i, sd->nalloc, "i", sopt->debug);
01401
01402 if ( ! rv && ! (sopt->skip_cols & V2S_SKIP_J) )
01403 rv = alloc_ints(&sd->j, sd->nalloc, "j", sopt->debug);
01404
01405 if ( ! rv && ! (sopt->skip_cols & V2S_SKIP_K) )
01406 rv = alloc_ints(&sd->k, sd->nalloc, "k", sopt->debug);
01407
01408 if ( ! rv && ! (sopt->skip_cols & V2S_SKIP_NVALS) )
01409 rv = alloc_ints(&sd->nvals, sd->nalloc, "nvals", sopt->debug);
01410
01411 if ( ! rv && ! (sopt->skip_cols & V2S_SKIP_VALS) )
01412 rv = alloc_vals_list(&sd->vals, sd->nalloc, sd->max_vals, sopt->debug);
01413 else if ( ! rv )
01414 rv = alloc_vals_list(&sd->vals, sd->nalloc, 1, sopt->debug);
01415
01416 if ( rv )
01417 {
01418 fprintf(stderr,"** failed v2s_results allocation\n");
01419 free_v2s_results(sd);
01420 sd = NULL;
01421 }
01422
01423 RETURN(sd);
01424 }
01425
01426
01427
01428
01429
01430
01431 static int alloc_vals_list(float *** ptr, int length, int width, int debug)
01432 {
01433 int c;
01434
01435 ENTRY("alloc_vals_list");
01436
01437 *ptr = (float **)malloc(width * sizeof(float *));
01438 if ( !*ptr )
01439 fprintf(stderr,"** avl: failed to alloc %d floats pointers\n", width);
01440
01441 for ( c = 0; c < width; c++ )
01442 {
01443 (*ptr)[c] = (float *)malloc(length * sizeof(float));
01444 if ( (*ptr)[c] == NULL )
01445 fprintf(stderr,"** avl: failed to alloc %d floats (# %d of %d)\n",
01446 length, c, width);
01447 }
01448
01449 if ( debug > 1 )
01450 fprintf(stderr,"-d alloc'd %d x %d floats for surf data\n",
01451 width, length);
01452
01453 RETURN(0);
01454 }
01455
01456
01457
01458
01459
01460
01461 static int realloc_vals_list(float ** ptr, int length, int width, int debug)
01462 {
01463 int c, count;
01464
01465 ENTRY("realloc_vals_list");
01466
01467 count = 0;
01468 for ( c = 0; c < width; c++ )
01469 {
01470 if ( ptr[c] )
01471 {
01472 ptr[c] = (float *)realloc(ptr[c], length * sizeof(float));
01473 if ( ptr[c] == NULL )
01474 fprintf(stderr,"** rvl: fail to realloc %d floats (%d of %d)\n",
01475 length, c, width);
01476 count++;
01477 }
01478 }
01479
01480 if ( debug > 1 )
01481 fprintf(stderr,"-d realloc'd %d x %d floats for surf data\n",
01482 count, length);
01483
01484 RETURN(0);
01485 }
01486
01487
01488
01489
01490
01491
01492 static int alloc_ints( int ** ptr, int length, char * dstr, int debug )
01493 {
01494 ENTRY("alloc_ints");
01495
01496 *ptr = (int *)malloc(length * sizeof(int));
01497 if ( ! *ptr )
01498 {
01499 fprintf(stderr,"** ai: failed to alloc %d ints for '%s'\n",length,dstr);
01500 RETURN(1);
01501 }
01502 if ( debug > 1 )
01503 fprintf(stderr,"-d ai: alloc'd %d ints for '%s'\n", length, dstr);
01504
01505 RETURN(0);
01506 }
01507
01508
01509
01510
01511
01512
01513 static int realloc_ints( int ** ptr, int length, char * dstr, int debug )
01514 {
01515 ENTRY("realloc_ints");
01516
01517 *ptr = (int *)realloc(*ptr, length * sizeof(int));
01518 if ( ! *ptr )
01519 {
01520 fprintf(stderr,"** ri: failed to alloc %d ints for '%s'\n",length,dstr);
01521 RETURN(1);
01522 }
01523 if ( debug > 1 )
01524 fprintf(stderr,"-d ri: alloc'd %d ints for '%s'\n", length, dstr);
01525
01526 RETURN(0);
01527 }
01528
01529
01530
01531
01532
01533
01534 int disp_v2s_param_t ( char * info, v2s_param_t * p )
01535 {
01536 ENTRY("disp_v2s_param_t");
01537
01538 if ( info )
01539 fputs( info, stderr );
01540
01541 if ( p == NULL )
01542 {
01543 fprintf(stderr, "disp_v2s_param_t: p == NULL\n" );
01544 RETURN(-1);
01545 }
01546
01547 fprintf(stderr,
01548 "v2s_param_t struct at %p :\n"
01549 " gpar : vcheck = %p : %s\n"
01550 " cmask = %p\n"
01551 " nvox, over_steps = %d, %d\n"
01552 " nsurf = %d\n"
01553 , p,
01554 p->gpar, ISVALID_DSET(p->gpar) ? "valid" : "invalid",
01555 p->cmask, p->nvox, p->over_steps, p->nsurf
01556 );
01557
01558 RETURN(0);
01559 }
01560
01561
01562
01563
01564
01565
01566 int disp_v2s_opts_t ( char * info, v2s_opts_t * sopt )
01567 {
01568 ENTRY("disp_v2s_opts_t");
01569
01570 if ( info )
01571 fputs( info, stderr );
01572
01573 if ( sopt == NULL )
01574 {
01575 fprintf(stderr, "disp_v2s_opts_t: sopt == NULL\n");
01576 RETURN(-1);
01577 }
01578
01579 fprintf(stderr,
01580 "v2s_opts_t struct at %p :\n"
01581 " map, gp_index = %d, %d\n"
01582 " debug, dnode = %d, %d\n"
01583 " no_head, skip_cols = %d, %d\n"
01584 " first_node, last_node = %d, %d\n"
01585 " use_norms, norm_len = %d, %f\n"
01586 " norm_dir = %d\n"
01587 " f_index, f_steps = %d, %d\n"
01588 " f_p1_fr, f_pn_fr = %f, %f\n"
01589 " f_p1_mm, f_pn_mm = %f, %f\n"
01590 " outfile_1D = %s\n"
01591 " outfile_niml = %s\n"
01592 , sopt,
01593 sopt->map, sopt->gp_index, sopt->debug, sopt->dnode,
01594 sopt->no_head, sopt->skip_cols,
01595 sopt->first_node, sopt->last_node,
01596 sopt->use_norms, sopt->norm_len, sopt->norm_dir,
01597 sopt->f_index, sopt->f_steps,
01598 sopt->f_p1_fr, sopt->f_pn_fr, sopt->f_p1_mm, sopt->f_pn_mm,
01599 CHECK_NULL_STR(sopt->outfile_1D),
01600 CHECK_NULL_STR(sopt->outfile_niml)
01601 );
01602
01603 RETURN(0);
01604 }
01605
01606
01607
01608
01609
01610
01611 static double magnitude_f( float * p, int length )
01612 {
01613 int c;
01614 double sums = 0.0;
01615
01616 for ( c = 0; c < length; c++, p++ )
01617 sums += (*p) * (*p);
01618
01619 return(sqrt(sums));
01620 }
01621
01622
01623
01624
01625
01626
01627 static float dist_f3mm( THD_fvec3 * p1, THD_fvec3 * p2 )
01628 {
01629 double d0, d1, d2;
01630
01631 if ( p1 == NULL || p2 == NULL )
01632 {
01633 fprintf( stderr, "** dist_f3mm: invalid params (%p,%p)\n", p1, p2 );
01634 return(0.0);
01635 }
01636
01637 d0 = p1->xyz[0] - p2->xyz[0];
01638 d1 = p1->xyz[1] - p2->xyz[1];
01639 d2 = p1->xyz[2] - p2->xyz[2];
01640
01641 return(sqrt(d0*d0 + d1*d1 + d2*d2));
01642 }
01643
01644
01645
01646
01647
01648
01649 static int disp_range_3dmm ( char * info, range_3dmm * dp )
01650 {
01651 ENTRY("disp_range_3dmm");
01652
01653 if ( info )
01654 fputs( info, stderr );
01655
01656 if ( dp == NULL )
01657 {
01658 fprintf(stderr, "disp_range_3dmm: dp == NULL\n");
01659 RETURN(-1);
01660 }
01661
01662 fprintf(stderr,
01663 "range_3dmm struct at %p :\n"
01664 " dset = %p : %s\n"
01665 " p1 = (%f, %f, %f)\n"
01666 " pn = (%f, %f, %f)\n"
01667 " dset_min = (%f, %f, %f)\n"
01668 " dset_max = (%f, %f, %f)\n"
01669 " oob_check, debug = (%d, %d)\n",
01670 dp, dp->dset, ISVALID_DSET(dp->dset) ? "valid" : "invalid",
01671 dp->p1.xyz[0], dp->p1.xyz[1], dp->p1.xyz[2],
01672 dp->pn.xyz[0], dp->pn.xyz[1], dp->pn.xyz[2],
01673 dp->dset_min.xyz[0], dp->dset_min.xyz[1], dp->dset_min.xyz[2],
01674 dp->dset_max.xyz[0], dp->dset_max.xyz[1], dp->dset_max.xyz[2],
01675 dp->oob_check, dp->debug
01676 );
01677
01678 RETURN(0);
01679 }
01680
01681
01682
01683
01684
01685
01686 static int disp_range_3dmm_res ( char * info, range_3dmm_res * dp )
01687 {
01688 ENTRY("disp_range_3dmm_res");
01689
01690 if ( info )
01691 fputs( info, stderr );
01692
01693 if ( dp == NULL )
01694 {
01695 fprintf(stderr, "disp_range_3dmm: dp == NULL\n");
01696 RETURN(-1);
01697 }
01698
01699 fprintf(stderr,
01700 "range_3dmm_res struct at %p :\n"
01701 " ims.imarr = %p\n"
01702 " ims.num, ims.nall = %d, %d\n"
01703 " repeats, masked = %d, %d\n"
01704 " oob, ifirst = %d, %d\n"
01705 " i3first[0,1,2] = %d, %d, %d\n"
01706 " i3arr = %p\n"
01707 , dp,
01708 dp->ims.imarr, dp->ims.num, dp->ims.nall,
01709 dp->repeats, dp->masked, dp->oob, dp->ifirst,
01710 dp->i3first.ijk[0], dp->i3first.ijk[1], dp->i3first.ijk[2],
01711 dp->i3arr );
01712
01713 if ( dp->i3arr )
01714 fprintf(stderr,
01715 " i3arr[0].ijk = %d, %d, %d\n",
01716 dp->i3arr[0].ijk[0], dp->i3arr[0].ijk[1], dp->i3arr[0].ijk[2] );
01717
01718 RETURN(0);
01719 }
01720
01721
01722
01723
01724
01725
01726 int disp_mri_imarr ( char * info, MRI_IMARR * dp )
01727 {
01728 float * fp;
01729 int cr, cc;
01730
01731 ENTRY("disp_mri_imarr");
01732
01733 if ( info )
01734 fputs( info, stderr );
01735
01736 if ( dp == NULL )
01737 {
01738 fprintf(stderr, "disp_mri_imarr: dp == NULL\n");
01739 RETURN(-1);
01740 }
01741
01742 fprintf(stderr,
01743 "mri_imarr struct at %p :\n"
01744 " num, nall = %d, %d\n",
01745 dp, dp->num, dp->nall );
01746
01747 for ( cr = 0; cr < dp->num; cr++ )
01748 {
01749 fp = MRI_FLOAT_PTR(dp->imarr[cr]);
01750 fprintf(stderr, " %3d: ", cr);
01751 for ( cc = 0; cc < dp->imarr[cr]->nx; cc++, fp++ )
01752 fprintf(stderr, "%f ", *fp );
01753 fputc( '\n', stderr );
01754 }
01755
01756 RETURN(0);
01757 }
01758
01759
01760
01761
01762
01763
01764 static int disp_surf_vals( char * mesg, v2s_results * sd, int node )
01765 {
01766 int index, c;
01767
01768 ENTRY("disp_surf_vals");
01769
01770 fprintf(stderr, "-------------------------------------------------\n");
01771 if ( mesg ) fputs( mesg, stderr );
01772 if ( sd->nused < 1 )
01773 {
01774 fprintf(stderr,"** no surf nodes defined\n");
01775 RETURN(-1);
01776 }
01777
01778 index = (node >= 0) ? node : sd->nused - 1;
01779
01780 fprintf(stderr, "v2s_results vals for sd_index %d, node %d :\n"
01781 " volind, (i, j, k) = %d, (%d, %d, %d)\n"
01782 " nvals: values... = %d: ", index,
01783 sd->nodes ? sd->nodes[index] : 0,
01784 sd->volind ? sd->volind[index] : 0,
01785 sd->i ? sd->i[index] : 0,
01786 sd->j ? sd->j[index] : 0,
01787 sd->k ? sd->k[index] : 0,
01788 sd->nvals ? sd->nvals[index] : 0 );
01789
01790 for ( c = 0; c < sd->max_vals; c++ )
01791 fprintf(stderr,"%s ", MV_format_fval(sd->vals[c][index]));
01792 fputc('\n', stderr);
01793
01794 RETURN(0);
01795 }
01796
01797
01798
01799
01800
01801
01802 int disp_v2s_results( char * mesg, v2s_results * d )
01803 {
01804 ENTRY("disp_v2s_results");
01805
01806 if ( mesg ) fputs( mesg, stderr );
01807
01808 fprintf(stderr, "v2s_results @ %p\n"
01809 " nalloc, nused = %d, %d\n"
01810 " max_vals, memory = %d, %d\n"
01811 " nodes, volind = %p, %p\n"
01812 " i, j, k = %p, %p, %p\n"
01813 " nvals, vals = %p, %p\n",
01814 d, d->nalloc, d->nused, d->max_vals, d->memory,
01815 d->nodes, d->volind, d->i, d->j, d->k, d->nvals, d->vals);
01816
01817 RETURN(0);
01818 }
01819
01820
01821
01822
01823
01824
01825 int disp_v2s_plugin_opts( char * mesg, v2s_plugin_opts * d )
01826 {
01827 ENTRY("disp_v2s_plugin_opts");
01828
01829 if ( mesg ) fputs( mesg, stderr );
01830
01831 fprintf(stderr, "v2s_plugin_opts @ %p\n"
01832 " ready = %d\n"
01833 " use0, use1 = %d, %d\n"
01834 " s0A, s0B = %d, %d\n"
01835 " s1A, s1B = %d, %d\n",
01836 d, d->ready, d->use0, d->use1,
01837 d->s0A, d->s0B, d->s1A, d->s1B );
01838 RETURN(0);
01839 }
01840
01841
01842
01843
01844
01845
01846 int free_v2s_results( v2s_results * sd )
01847 {
01848 int c;
01849
01850 ENTRY("free_v2s_results");
01851
01852 if( ! sd ) RETURN(0);
01853
01854 if (sd->nodes) { free(sd->nodes); sd->nodes = NULL; }
01855 if (sd->volind) { free(sd->volind); sd->volind = NULL; }
01856 if (sd->i) { free(sd->i); sd->i = NULL; }
01857 if (sd->j) { free(sd->j); sd->j = NULL; }
01858 if (sd->k) { free(sd->k); sd->k = NULL; }
01859 if (sd->nvals) { free(sd->nvals); sd->nvals = NULL; }
01860
01861 if (sd->vals)
01862 {
01863 for ( c = 0; c < sd->max_vals; c++ )
01864 if ( sd->vals[c] ) { free(sd->vals[c]); sd->vals[c] = NULL; }
01865
01866 free(sd->vals);
01867 sd->vals = NULL;
01868 }
01869
01870 free(sd);
01871
01872 RETURN(0);
01873 }
01874
01875
01876
01877
01878
01879
01880 static int validate_v2s_inputs ( v2s_opts_t * sopt, v2s_param_t * p )
01881 {
01882 int c;
01883
01884 ENTRY("validate_v2s_inputs");
01885
01886 if ( !sopt || !p )
01887 {
01888 fprintf(stderr,"** validate_v2s_inputs: bad params (%p,%p)\n", sopt, p);
01889 RETURN(-1);
01890 }
01891
01892
01893 if ( sopt->map <= E_SMAP_INVALID || sopt->map >= E_SMAP_FINAL )
01894 {
01895 fprintf(stderr,"** invalid mapping index %d\n", sopt->map);
01896 RETURN(1);
01897 }
01898 else if ( v2s_map_type(gv2s_map_names[sopt->map]) != sopt->map )
01899 {
01900 fprintf(stderr,"** mapping index failure for %d\n", sopt->map);
01901 RETURN(1);
01902 }
01903
01904 if ( sopt->f_steps <= 0 || sopt->f_steps >= V2S_STEPS_TOOOOO_BIG )
01905 {
01906 fprintf(stderr,"** invalid f_steps = %d\n", sopt->f_steps);
01907 RETURN(1);
01908 }
01909
01910
01911 if ( !p->gpar ) {fprintf(stderr,"** vv2si: no dset?\n"); RETURN(2);}
01912
01913 if ( p->nvox != DSET_NVOX(p->gpar) )
01914 {
01915 fprintf(stderr,"** invalid voxel count (%d) for grid_parent\n",p->nvox);
01916 RETURN(2);
01917 }
01918
01919 if ( sopt->gp_index >= DSET_NVALS(p->gpar) )
01920 {
01921 fprintf(stderr,"** gp_index (%d) > max grid_parent index (%d)\n",
01922 sopt->gp_index, DSET_NVALS(p->gpar) - 1);
01923 RETURN(2);
01924 }
01925
01926 if ( p->nsurf < 1 || p->nsurf > 2 )
01927 {
01928 fprintf(stderr,"** invalid: nsurf = %d must be in [%d,%d]\n",
01929 p->nsurf, 1, 2);
01930 RETURN(2);
01931 }
01932
01933
01934 for ( c = 0; c < p->nsurf; c++ )
01935 if ( check_SUMA_surface(p->surf + c) )
01936 RETURN(3);
01937
01938
01939 if ( p->nsurf > 1 )
01940 {
01941 if ( p->surf[0].num_ixyz != p->surf[1].num_ixyz )
01942 {
01943 fprintf(stderr,"** invalid surfaces, different # nodes (%d,%d)\n",
01944 p->surf[0].num_ixyz, p->surf[1].num_ixyz);
01945 RETURN(4);
01946 }
01947 }
01948 else if ( sopt->use_norms && !p->surf[0].norm )
01949 {
01950 fprintf(stderr,"** missing surface normals for surface '%s'\n",
01951 CHECK_EMPTY_STR(p->surf[0].label));
01952 RETURN(4);
01953 }
01954
01955 RETURN(0);
01956 }
01957
01958
01959 static int check_SUMA_surface( SUMA_surface * s )
01960 {
01961 int rv = 0;
01962 ENTRY("is_valid_SUMA_surface");
01963
01964 if ( !s ) { fprintf(stderr,"** ivSs: no surface\n"); RETURN(0); }
01965
01966 if ( s->type != SUMA_SURFACE_TYPE )
01967 {
01968 fprintf(stderr,"** surf '%s': invalid type %d\n",
01969 CHECK_EMPTY_STR(s->label), s->type);
01970 rv++;
01971 }
01972
01973 if ( s->num_ixyz < 0 || s->nall_ixyz < 0 || s->num_ixyz < s->nall_ixyz )
01974 {
01975 fprintf(stderr,"** surf '%s': invalid num_ixyz = %d, nall_ixyz = %d\n",
01976 CHECK_EMPTY_STR(s->label), s->num_ixyz, s->nall_ixyz);
01977 rv++;
01978 }
01979
01980 if ( s->seq == 0 || s->sorted == 0 || s->seqbase != 0 )
01981 {
01982 fprintf(stderr,"** surf '%s': invalid seq, sort or base (%d, %d, %d)\n",
01983 CHECK_EMPTY_STR(s->label), s->seq, s->sorted, s->seqbase);
01984 rv++;
01985 }
01986
01987 if ( !s->ixyz )
01988 {
01989 fprintf(stderr,"** surf '%s': invalid, missing nodes\n",
01990 CHECK_EMPTY_STR(s->label));
01991 rv++;
01992 }
01993
01994 RETURN(rv);
01995 }
01996
01997
01998
01999
02000
02001
02002
02003
02004
02005 int v2s_vals_over_steps( int map )
02006 {
02007 if ( map == E_SMAP_SEG_VALS )
02008 return 1;
02009
02010 return 0;
02011 }
02012
02013
02014
02015
02016
02017
02018
02019
02020
02021 int v2s_map_type ( char * map_str )
02022 {
02023 v2s_map_nums map;
02024
02025 ENTRY("v2s_map_type");
02026
02027 if ( map_str == NULL )
02028 {
02029 fprintf( stderr, "** v2s_map_type: missing map_str parameter\n" );
02030 RETURN((int)E_SMAP_INVALID);
02031 }
02032
02033 if ( sizeof(gv2s_map_names) / sizeof(char *) != (int)E_SMAP_FINAL )
02034 {
02035 fprintf( stderr, "** error: gv2s_map_names/v2s_map_num mis-match\n");
02036 RETURN((int)E_SMAP_INVALID);
02037 }
02038
02039
02040 if ( !strcmp( map_str, gv2s_map_names[E_SMAP_COUNT] ) )
02041 RETURN((int)E_SMAP_INVALID);
02042
02043 for ( map = E_SMAP_INVALID; map < E_SMAP_FINAL; map++ )
02044 if ( !strcmp( map_str, gv2s_map_names[map] ) )
02045 RETURN((int)map);
02046
02047 RETURN((int)E_SMAP_INVALID);
02048 }
02049
02050
02051
02052
02053
02054
02055
02056
02057 int thd_mask_from_brick(THD_3dim_dataset * dset, int volume, float thresh,
02058 byte ** mask, int absolute)
02059 {
02060 float factor;
02061 byte * tmask;
02062 int nvox, type, c, size = 0;
02063
02064 ENTRY("thd_mask_from_brick");
02065
02066 if ( mask ) *mask = NULL;
02067
02068 if ( !ISVALID_DSET(dset) || ! mask || volume < 0 )
02069 RETURN(-1);
02070
02071 if ( volume >= DSET_NVALS(dset) )
02072 {
02073 fprintf(stderr,"** tmfb: sub-brick %d out-of-range\n", volume);
02074 RETURN(-1);
02075 }
02076
02077 nvox = DSET_NVOX(dset);
02078 type = DSET_BRICK_TYPE(dset, volume);
02079
02080 if ( type != MRI_byte && type != MRI_short &&
02081 type != MRI_int && type != MRI_float )
02082 {
02083 fprintf(stderr,"** tmfb: invalid dataset type %s, sorry...\n",
02084 MRI_type_name[type]);
02085 RETURN(-1);
02086 }
02087
02088 tmask = (byte *)calloc(nvox, sizeof(byte));
02089 if ( ! tmask )
02090 {
02091 fprintf(stderr,"** tmfb: failed to allocate mask of %d bytes\n", nvox);
02092 RETURN(-1);
02093 }
02094
02095 factor = DSET_BRICK_FACTOR(dset, volume);
02096
02097
02098 if ( factor != 0.0 ) thresh /= factor;
02099
02100 switch( DSET_BRICK_TYPE(dset, volume) )
02101 {
02102 case MRI_byte:
02103 {
02104 byte * dp = DSET_ARRAY(dset, volume);
02105 byte thr = BYTEIZE(thresh + 0.99999);
02106 for ( c = 0; c < nvox; c++ )
02107 if ( dp[c] != 0 && ( dp[c] >= thr ) )
02108 {
02109 size++;
02110 tmask[c] = 1;
02111 }
02112 }
02113 break;
02114
02115 case MRI_short:
02116 {
02117 short * dp = DSET_ARRAY(dset, volume);
02118 short thr = SHORTIZE(thresh + 0.99999);
02119 for ( c = 0; c < nvox; c++, dp++ )
02120 if ( *dp != 0 && ( *dp >= thr || (absolute && *dp <= -thr) ) )
02121 {
02122 size++;
02123 tmask[c] = 1;
02124 }
02125 }
02126 break;
02127
02128 case MRI_int:
02129 {
02130 int * dp = DSET_ARRAY(dset, volume);
02131 int thr = (int)(thresh + 0.99999);
02132 for ( c = 0; c < nvox; c++, dp++ )
02133 if ( *dp != 0 && ( *dp >= thr || (absolute && *dp <= -thr) ) )
02134 {
02135 size++;
02136 tmask[c] = 1;
02137 }
02138 }
02139 break;
02140
02141 case MRI_float:
02142 {
02143 float * dp = DSET_ARRAY(dset, volume);
02144 for ( c = 0; c < nvox; c++, dp++ )
02145 if (*dp != 0 && (*dp >= thresh || (absolute && *dp <= -thresh)))
02146 {
02147 size++;
02148 tmask[c] = 1;
02149 }
02150 }
02151 break;
02152
02153 default:
02154 {
02155 fprintf(stderr,"** tmfb: invalid dataset type, sorry...\n");
02156 free(tmask);
02157 }
02158 break;
02159 }
02160
02161 *mask = tmask;
02162
02163 RETURN(size);
02164 }
02165
02166
02167
02168
02169
02170
02171
02172
02173 int v2s_is_good_map( int map, int from_afni )
02174 {
02175 ENTRY("v2s_good_map_index");
02176
02177 if ( map <= E_SMAP_INVALID || map >= E_SMAP_FINAL )
02178 RETURN(0);
02179
02180
02181 if ( map == E_SMAP_MASK2 || map == E_SMAP_COUNT )
02182 RETURN(0);
02183
02184 if ( from_afni && map == E_SMAP_SEG_VALS )
02185 RETURN(0);
02186
02187 RETURN(1);
02188 }
02189
02190
02191
02192
02193
02194
02195
02196
02197 static int fill_sopt_default(v2s_opts_t * sopt, int nsurf )
02198 {
02199
02200 ENTRY("fill_sopt_default");
02201
02202 if ( !sopt || nsurf < 1 || nsurf > 2 )
02203 {
02204 fprintf(stderr,"** fill_sopt_default: bad params (%p,%d)\n",sopt,nsurf);
02205 RETURN(1);
02206 }
02207
02208
02209 memset(sopt, 0, sizeof(*sopt));
02210
02211 if ( nsurf == 2 )
02212 sopt->map = E_SMAP_MIDPT;
02213 else
02214 sopt->map = E_SMAP_MASK;
02215
02216 sopt->gp_index = -1;
02217 sopt->no_head = 1;
02218 sopt->skip_cols = V2S_SKIP_ALL ^ V2S_SKIP_NODES;
02219 sopt->f_steps = 1;
02220 sopt->outfile_1D = NULL;
02221 sopt->outfile_niml = NULL;
02222
02223 RETURN(0);
02224 }
02225
02226
02227
02228
02229
02230
02231 int v2s_write_outfile_1D ( v2s_opts_t * sopt, v2s_results * sd, char * label )
02232 {
02233 FILE * fp;
02234 int c, c2;
02235 ENTRY("v2s_write_outfile_1D");
02236
02237 fp = fopen( sopt->outfile_1D, "w" );
02238 if ( fp == NULL )
02239 {
02240 fprintf( stderr, "** failure to open '%s' for writing\n",
02241 sopt->outfile_1D );
02242 RETURN(-1);
02243 }
02244
02245 if ( ! sopt->no_head )
02246 print_header(fp, label, gv2s_map_names[sopt->map], sd);
02247
02248 for ( c = 0; c < sd->nused; c++ )
02249 {
02250
02251 fputc(' ', fp);
02252 if ( sd->nodes ) fprintf(fp, " %8d", sd->nodes[c]);
02253 if ( sd->volind ) fprintf(fp, " %8d ", sd->volind[c]);
02254 if ( sd->i ) fprintf(fp, " %3d", sd->i[c]);
02255 if ( sd->j ) fprintf(fp, " %3d", sd->j[c]);
02256 if ( sd->k ) fprintf(fp, " %3d", sd->k[c]);
02257 if ( sd->nvals ) fprintf(fp, " %3d", sd->nvals[c]);
02258
02259 for ( c2 = 0; c2 < sd->max_vals; c2++ )
02260 fprintf(fp, " %10s", MV_format_fval(sd->vals[c2][c]));
02261 fputc('\n', fp);
02262 }
02263
02264 fclose(fp);
02265
02266 RETURN(0);
02267 }
02268
02269
02270
02271
02272
02273
02274
02275 int v2s_write_outfile_niml ( v2s_opts_t * sopt, v2s_results * sd, int free_vals){
02276 static char v2s_name[] = "3dVol2Surf_dataset";
02277 NI_element * nel = NULL;
02278 NI_stream ns;
02279 char * ni_name;
02280 int c;
02281
02282 ENTRY("v2s_write_outfile_niml");
02283
02284 if ( !sopt->outfile_niml ) RETURN(0);
02285
02286 nel = NI_new_data_element( v2s_name, sd->nused );
02287 if ( !nel )
02288 {
02289 fprintf(stderr,"** file NI_new_data_element, n = '%s', len = %d\n",
02290 v2s_name, sd->nused);
02291 RETURN(1);
02292 }
02293
02294 ni_name = (char *)calloc(strlen(sopt->outfile_niml)+6, sizeof(char));
02295 if ( !ni_name ) { fprintf(stderr,"** ni_name failed\n"); RETURN(1); }
02296 sprintf(ni_name, "file:%s", sopt->outfile_niml);
02297
02298 ns = NI_stream_open(ni_name, "w");
02299
02300 NI_add_column(nel,NI_INT,sd->nodes);
02301
02302 for ( c = 0; c < sd->max_vals; c++ )
02303 {
02304 NI_add_column(nel, NI_FLOAT, sd->vals[c]);
02305 if ( free_vals ) { free(sd->vals[c]); sd->vals[c] = NULL; }
02306 }
02307 if ( free_vals ) { free(sd->vals); sd->vals = NULL; }
02308
02309 if ( NI_write_element(ns, nel, NI_BINARY_MODE) < 0 )
02310 {
02311 fprintf(stderr,"** NI_write_element failed for: '%s'\n", ni_name);
02312 RETURN(1);
02313 }
02314
02315 NI_free_element( nel );
02316 NI_stream_close( ns );
02317 free(ni_name);
02318
02319 RETURN(0);
02320 }
02321
02322
02323
02324
02325
02326
02327 static int print_header(FILE * outfp, char * surf, char * map, v2s_results * sd)
02328 {
02329 int val;
02330
02331 ENTRY("print_header");
02332
02333 fprintf( outfp, "# --------------------------------------------------\n" );
02334 fprintf( outfp, "# surface '%s', '%s' :\n", surf, map );
02335 fprintf( outfp, "#\n" );
02336
02337
02338
02339
02340
02341
02342
02343 fputc( '#', outfp );
02344
02345 if ( sd->nodes ) fprintf(outfp, " node ");
02346 if ( sd->volind ) fprintf(outfp, " 1dindex ");
02347 if ( sd->i ) fprintf(outfp, " i ");
02348 if ( sd->j ) fprintf(outfp, " j ");
02349 if ( sd->k ) fprintf(outfp, " k ");
02350 if ( sd->nvals ) fprintf(outfp, " vals");
02351
02352 for ( val = 0; val < sd->max_vals; val++ )
02353 fprintf( outfp, " v%-2d ", val );
02354 fputc( '\n', outfp );
02355
02356 fputc( '#', outfp );
02357
02358 if ( sd->nodes ) fprintf(outfp, " ------");
02359 if ( sd->volind ) fprintf(outfp, " ------- ");
02360 if ( sd->i ) fprintf(outfp, " ---");
02361 if ( sd->j ) fprintf(outfp, " ---");
02362 if ( sd->k ) fprintf(outfp, " ---");
02363 if ( sd->nvals ) fprintf(outfp, " ----");
02364
02365 fputs( " ", outfp );
02366 for ( val = 0; val < sd->max_vals; val++ )
02367 fprintf( outfp, " -------- " );
02368 fputc( '\n', outfp );
02369
02370 RETURN(0);
02371 }
02372
02373
