00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 #include "afni.h"
00020 #include "plug_maxima.h"
00021
00022 #ifndef ALLOW_PLUGINS
00023 # error "Plugins not properly set up -- see machdep.h"
00024 #endif
00025
00026
00027
00028
00029
00030 char * MAXIMA_main( PLUGIN_interface * );
00031
00032 char grMessage[ R_MESSAGE_L ];
00033 static char * grStyle[] = { "Sort-n-Remove", "Weighted-Average" };
00034 static char * grSvals[] = { "1 (default)", "1 to N", "N to 1" };
00035 static char * grNY[] = { "No", "Yes" };
00036
00037 static char helpstring[] =
00038 "Maxima - used to locate extrema in a functional dataset.\n"
00039 "\n"
00040 "This plugin reads a functional dataset and locates any relative extrema\n"
00041 "(maximums or minimums, depending on the user option). A _relative_\n"
00042 "maximum is a point that is greater than all neighbors (not necessarily\n"
00043 "greater than all other values in the sub-brick). The output from this\n"
00044 "process can be text based (sent to the terminal window) and it can be a\n"
00045 "mask (integral) dataset, where the locations of the extrema are set.\n"
00046 "\n"
00047 "When writing a dataset, it is often useful to set a sphere around each\n"
00048 "extrema, not to just set individual voxels. This makes viewing those\n"
00049 "locations much more reasonable. Also, if the 'Sphere Values' option is\n"
00050 "set to 'N to 1', the sphere around the most extreme voxel will get the\n"
00051 "value N, giving it the 'top' color in afni (and so on, down to 1).\n"
00052 "\n"
00053 "Notes : The only required option is the input dataset.\n"
00054 " Input datasets must be of type short.\n"
00055 " All distances are in voxel units.\n"
00056 "\n"
00057 " *** Options ***\n"
00058 "\n"
00059 "----- Input: -----\n"
00060 "\n"
00061 " Dataset - dataset to locate extrema within\n"
00062 "\n"
00063 " Sub-brick - sub-brick index of volume to locate extrema within\n"
00064 "\n"
00065 "----- Output Dset: -----\n"
00066 "\n"
00067 " Prefix - prefix for an output functional dataset\n"
00068 "\n"
00069 " This dataset may be viewed as a mask. Its set points\n"
00070 " will correspond to any selected extrema. The \"Output\n"
00071 " Size\" option will change those points to spheres.\n"
00072 "\n"
00073 " Values - values to fill output spheres with\n"
00074 "\n"
00075 " 1 (default) - fill all spheres with the value 1\n"
00076 " 1 to N - fill sphere of most extreme point with 1,\n"
00077 " the next with 2, and so on\n"
00078 " N to 1 - fill sphere of most extreme point with N,\n"
00079 " the next with N-1, etc., perhaps matching\n"
00080 " Overlay coloring better\n"
00081 "\n"
00082 "----- Threshold: -----\n"
00083 "\n"
00084 " Cutoff - provides a cutoff value for extrema\n"
00085 "\n"
00086 " Extrema not meeting this cutoff will not be considered.\n"
00087 " Note that if the \"Neg Extrema\" option is applied, the\n"
00088 " user will generally want a negative threshold.\n"
00089 "\n"
00090 "----- Separation: -----\n"
00091 "\n"
00092 " Distance - minimum acceptable distance between extrema\n"
00093 "\n"
00094 " Less significant extrema which are too close to more\n"
00095 " significant extrema will be discounted in some way.\n"
00096 " This option works with the \"Neighbor Style\" option.\n"
00097 "\n"
00098 " Note that the distance is in voxels, not mm.\n"
00099 "\n"
00100 "----- Output Size: -----\n"
00101 "\n"
00102 " Radius - output mask radius around set points\n"
00103 "\n"
00104 " If the user wants the output BRIK to consist of spheres\n"
00105 " centered around extrema points, this option can be used\n"
00106 " to set the radius for those spheres.\n"
00107 "\n"
00108 "----- Neighbor: -----\n"
00109 "\n"
00110 " Style - technique for handling clustered extrema\n"
00111 "\n"
00112 " If extrema are not as far apart as is specified by the\n"
00113 " \"Separation\" option, the Neighbor (Style) option\n"
00114 " specifies how to handle those points.\n"
00115 "\n"
00116 " Sort-n-Remove : The extrema are sorted by magnitude.\n"
00117 " For each such extrema, if the extrema is still set\n"
00118 " (not previously removed), all (less significant)\n"
00119 " neighbors within the Separation radius are removed.\n"
00120 "\n"
00121 " Weighted-Average : Again, traverse the sorted list of\n"
00122 " extrema. Set the current extrema to the center of\n"
00123 " mass of all extrema within the Separation radius.\n"
00124 "\n"
00125 "----- Params: -----\n"
00126 "\n"
00127 " Neg Extr - search for negative extrema (minima)\n"
00128 "\n"
00129 " This will search for the minima of the dataset.\n"
00130 " Note that a negative threshold may be desired.\n"
00131 "\n"
00132 " True Max - do not consider extrema with an equal valued neighbor\n"
00133 "\n"
00134 " As a default points may be considered extrema, even if\n"
00135 " they have a neighbor of the same value. Setting this\n"
00136 " option requires extrema to be stricktly greater than\n"
00137 " any of their neighbors.\n"
00138 " Note that with this set, extrema locations that have\n"
00139 " more than one voxel at the maximal value are ignored.\n"
00140 "\n"
00141 "----- Output Text: -----\n"
00142 "\n"
00143 " No Text - do not display the extrma points as text\n"
00144 "\n"
00145 " Debug - how much extra information to output to the terminal\n"
00146 "\n"
00147 "----- Output Coords: -----\n"
00148 "\n"
00149 " Dicom - Yes/No: whether to display output in Dicom orientation\n"
00150 "\n"
00151 " In Dicom orientation, the output is RAI, meaning right,\n"
00152 " anterior and inferior are negative coordinates, and they\n"
00153 " are printed in that order (RL, then AP, then IS).\n"
00154 "\n"
00155 " If this option is set to NO, the dataset orientation is\n"
00156 " used, whichever of the 48 it happens to be.\n"
00157 "\n"
00158 " Note that in either case, the output orientation is\n"
00159 " printed above the results, in the terminal window.\n"
00160 "\n"
00161 "Author: (the hopefully-no-longer-lamented) R Reynolds\n"
00162 "\n"
00163 "History:\n"
00164 "\n"
00165 " 20 Feb 2004 [rickr]\n"
00166 " - added ENTRY/RETURN macros\n"
00167 " - do not process last plane\n"
00168 " - allow any anat/func type datasets of type short\n"
00169 " - added sub-brick selector\n"
00170 " - note that dist and radius are in voxels\n"
00171 " - output coordinates in RAI mm format\n"
00172 "\n"
00173 " 01 Nov 2004 [rickr]\n"
00174 " - remove restrictions on threshold input\n"
00175 " - rearrange options, and add a Debug Level\n"
00176 " - increment style (should be in {1,2}, not {0,1}\n"
00177 " - add a little debug output, including show_point_list_s()\n"
00178 " - removed unused variables\n"
00179 " - true_max update in find_local_maxima()\n"
00180 " - added check for warp-on-demand failure\n"
00181 "\n"
00182 " 14 Feb 2005 [rickr]\n"
00183 " - added the 'Sphere Values' and 'Dicom Coords' options\n"
00184 "\n"
00185 " 07 Mar 2005 [rickr]\n"
00186 " - output appropriate coords via new THD_3dind_to_3dmm_no_wod()\n"
00187 " - added new debug output\n"
00188 " - changed default separation to 4 voxels\n"
00189 " - added gr_fac for printing data values in debug mode\n"
00190 ;
00191
00192
00193
00194
00195
00196
00197
00198
00199 DEFINE_PLUGIN_PROTOTYPE
00200
00201 PLUGIN_interface * PLUGIN_init( int ncall )
00202 {
00203 PLUGIN_interface * plint ;
00204
00205 if( ncall > 0 ) return NULL ;
00206
00207
00208
00209 plint = PLUTO_new_interface( "Maxima", "find extrema in a dataset",
00210 helpstring, PLUGIN_CALL_VIA_MENU , MAXIMA_main );
00211
00212 PLUTO_add_hint( plint, "find local maxima/minima" );
00213
00214 PLUTO_set_sequence( plint , "z:Reynolds" ) ;
00215
00216
00217
00218 PLUTO_add_option( plint, "Input" , "Input" , TRUE );
00219 PLUTO_add_hint( plint, "choose dataset for input" );
00220 PLUTO_add_dataset(plint, "Dataset" , ANAT_ALL_MASK , FUNC_ALL_MASK,
00221 DIMEN_ALL_MASK | BRICK_SHORT_MASK );
00222 PLUTO_add_number( plint , "Sub-brick" , 0,9999,0 , 0,1 ) ;
00223
00224
00225
00226
00227 PLUTO_add_option( plint, "Output Dset" , "prefix" , FALSE );
00228 PLUTO_add_hint( plint, "options for the creation of an output dataset");
00229 PLUTO_add_string( plint, "Prefix", 0 , NULL, 19 );
00230 PLUTO_add_hint( plint, "option: choose dataset prefix for output" );
00231 PLUTO_add_string( plint, "Sphere Values", 3 , grSvals, 0 );
00232 PLUTO_add_hint( plint, "option: choose value style for output spheres" );
00233
00234
00235
00236 PLUTO_add_option( plint, "Threshold" , "cutoff" , FALSE ) ;
00237 PLUTO_add_hint( plint, "option: choose a threshold for value at extrema" );
00238 PLUTO_add_number( plint, "Cutoff", 0, 0, 0, 1000, 1 );
00239
00240
00241
00242 PLUTO_add_option( plint, "Separation" , "min_dist" , FALSE ) ;
00243 PLUTO_add_hint( plint, "option: choose a minimum distance between extrema" );
00244 PLUTO_add_number( plint, "Distance(vox)", 0, 1000, 1, 40, 1 );
00245
00246
00247
00248 PLUTO_add_option( plint, "Output Size" , "out_rad" , FALSE ) ;
00249 PLUTO_add_hint( plint, "option: choose a spherical radius around extrema "
00250 "points in mask" );
00251 PLUTO_add_number( plint, "Radius(vox)", 0, 1000, 1, 50, 1 );
00252
00253
00254
00255 PLUTO_add_option( plint, "Neighbor" , "style" , FALSE ) ;
00256 PLUTO_add_hint( plint, "option: technique for neighbor removal" );
00257 PLUTO_add_string( plint, "Style", 2, grStyle, 0 );
00258
00259
00260
00261 PLUTO_add_option( plint, "Params" , "params" , FALSE ) ;
00262 PLUTO_add_hint( plint, "options: negative extrema and true max" );
00263 PLUTO_add_string( plint, "Neg Extrema", 2, grNY, 0 );
00264 PLUTO_add_hint( plint, "search for negative extrema, not positive" );
00265 PLUTO_add_string( plint, "True Max", 2, grNY, 0 );
00266 PLUTO_add_hint( plint, "exclude extrema with equal neighbors" );
00267
00268
00269
00270 PLUTO_add_option( plint, "Output Text" , "output" , FALSE ) ;
00271 PLUTO_add_hint( plint, "options: no output text, debug level" );
00272 PLUTO_add_string( plint, "No Text Out", 2, grNY, 0 );
00273 PLUTO_add_hint( plint, "do not output extrema as text (to terminal)" );
00274 PLUTO_add_number( plint, "Debug Level", 0, 4, 0, 0, 0 );
00275 PLUTO_add_hint( plint, "search for negative extrema, not positive" );
00276
00277
00278
00279 PLUTO_add_option( plint, "Output Coords" , "dicom_coords" , FALSE ) ;
00280 PLUTO_add_hint( plint, "option: output coordinates in Dicom format" );
00281 PLUTO_add_string( plint, "Dicom Coords", 2, grNY, 1 );
00282
00283 return plint;
00284 }
00285
00286
00287
00288
00289
00290
00291
00292
00293 char * MAXIMA_main( PLUGIN_interface * plint )
00294 {
00295 r_afni_s A;
00296 maxima_s M;
00297 char * ret_string = NULL;
00298
00299
00300 if ( ( ret_string = process_args( &A, &M, plint ) ) != NULL )
00301 return ret_string;
00302
00303 if ( ! process_data( &M ) )
00304 return "************************************\n"
00305 "MAXIMA_main: data processing failure\n"
00306 "************************************";
00307
00308 if ( ! write_results( &A, &M, plint ) )
00309 return "***********************************\n"
00310 "MAXIMA_main: result writing failure\n"
00311 "***********************************";
00312
00313 free_memory( &A, &M );
00314
00315 return NULL;
00316 }
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327 static char *
00328 process_args( r_afni_s * A, maxima_s * M, PLUGIN_interface * plint )
00329 {
00330 THD_3dim_dataset * dset;
00331 MCW_idcode * idc ;
00332 char * optag, * outfile = NULL, * str;
00333 float cutoff = 0.0, min_dist = 0.0, out_rad = 0.0;
00334 int negatives = 0, quiet = 0, true_max = 0, opcnt = 0;
00335 int val, debug = 0, style = MAX_SORT_N_REMOVE_STYLE, sb;
00336 int sval_style = 0, dicom_coords = 1;
00337
00338 ENTRY("process_args");
00339
00340
00341 if( plint == NULL )
00342 RETURN("----------------------\n"
00343 "arguments : NULL input\n"
00344 "----------------------");
00345
00346 if ( ! init_afni_s( A ) )
00347 RETURN( "------------------------\n"
00348 "arguments : init failure\n"
00349 "------------------------");
00350
00351 PLUTO_next_option( plint );
00352 idc = PLUTO_get_idcode( plint );
00353 dset = PLUTO_find_dset( idc );
00354
00355 if( dset == NULL )
00356 RETURN("-----------------------------\n"
00357 "arguments : bad input dataset\n"
00358 "-----------------------------");
00359
00360 sb = (int)PLUTO_get_number( plint );
00361 if ( sb >= DSET_NVALS(dset) || sb < 0 )
00362 RETURN("--------------------------\n"
00363 "arguments : bad sub-brick \n"
00364 "--------------------------");
00365 A->sub_brick = sb;
00366
00367 DSET_load( dset );
00368
00369 for ( optag = PLUTO_get_optiontag( plint );
00370 optag != NULL;
00371 optag = PLUTO_get_optiontag( plint )
00372 )
00373 {
00374 if ( ! strcmp( optag, "prefix" ) )
00375 {
00376 outfile = PLUTO_get_string( plint );
00377 if ( ! PLUTO_prefix_ok( outfile ) )
00378 RETURN( "-------------------------\n"
00379 "options : bad file prefix\n"
00380 "-------------------------");
00381 sval_style = PLUTO_string_index(PLUTO_get_string(plint),3,grSvals);
00382 }
00383 else if ( ! strcmp( optag, "cutoff" ) )
00384 {
00385 cutoff = PLUTO_get_number( plint );
00386 }
00387 else if ( ! strcmp( optag, "min_dist" ) )
00388 {
00389 if ( ( min_dist = PLUTO_get_number( plint ) ) < 0 )
00390 RETURN( "-----------------------------------------\n"
00391 "options : Separation must be non-negative\n"
00392 "-----------------------------------------");
00393 }
00394 else if ( ! strcmp( optag, "out_rad" ) )
00395 {
00396 if ( ( out_rad = PLUTO_get_number( plint ) ) < 0 )
00397 RETURN( "--------------------------------------------\n"
00398 "options : Output radius must be non-negative\n"
00399 "--------------------------------------------");
00400 }
00401 else if ( ! strcmp( optag, "params" ) )
00402 {
00403 str = PLUTO_get_string( plint );
00404 val = PLUTO_string_index(str, 2, grNY);
00405 if ( val > 0 ) negatives = 1;
00406 str = PLUTO_get_string( plint );
00407 val = PLUTO_string_index(str, 2, grNY);
00408 if ( val > 0 ) true_max = 1;
00409 }
00410 else if ( ! strcmp( optag, "style" ) )
00411 {
00412 if ( ( str = PLUTO_get_string( plint ) ) == NULL )
00413 RETURN( "-------------------------------\n"
00414 "options : missing style string?\n"
00415 "-------------------------------");
00416 if ((( style = PLUTO_string_index(str, 2, grStyle))
00417 < 0 ) || ( style >= MAX_MAX_STYLE ) )
00418 {
00419 sprintf( grMessage,
00420 "---------------------------\n"
00421 "options : bad style is %d\n"
00422 "---------------------------", style );
00423 RETURN(grMessage);
00424 }
00425 style++;
00426 }
00427 else if ( ! strcmp( optag, "output" ) )
00428 {
00429 str = PLUTO_get_string( plint );
00430 val = PLUTO_string_index(str, 2, grNY);
00431 if ( val > 0 ) quiet = 1;
00432 debug = PLUTO_get_number( plint );
00433
00434 }
00435 else if ( ! strcmp( optag, "dicom_coords" ) )
00436 {
00437 str = PLUTO_get_string( plint );
00438 val = PLUTO_string_index(str, 2, grNY);
00439 if ( val == 0 ) dicom_coords = 0;
00440 }
00441 else
00442 {
00443 sprintf( grMessage, "Error: pa_00\n"
00444 "Unexpected option #%d: '%s'", opcnt, optag );
00445 RETURN(grMessage);
00446 }
00447
00448 opcnt++;
00449 }
00450
00451 if ( ( out_rad > 0 ) && ( outfile == NULL ) )
00452 RETURN( "------------------------------------------------\n"
00453 "arguments : specify outfile to use output radius\n"
00454 "------------------------------------------------");
00455
00456 if ( ! r_set_afni_s_from_dset( A, dset ) )
00457 RETURN( "-------------------------------\n"
00458 "arguments : afni_s init failure\n"
00459 "-------------------------------");
00460
00461 if ( ! init_maxima_s( M, A, outfile ) )
00462 RETURN("----------------------------------\n"
00463 "MAXIMA_main: maxima_s init failure\n"
00464 "----------------------------------");
00465
00466
00467 M->sval_style = sval_style;
00468 M->cutoff = cutoff / A->factor[0];
00469 M->min_dist = min_dist;
00470 M->out_rad = out_rad;
00471
00472 M->negatives = negatives;
00473 M->ngbr_style = style;
00474 M->quiet = quiet;
00475 M->true_max = true_max;
00476 M->dicom_coords = dicom_coords;
00477 M->debug = debug;
00478
00479 gr_fac = A->factor[0];
00480
00481 if ( M->debug > 0 )
00482 {
00483 if ( M->debug > 3 ) show_maxima_s( "plugin values applied ", M );
00484 fprintf(stderr," using sub-brick %d, factor %f (1/%f)\n",
00485 A->sub_brick, A->factor[0], 1/A->factor[0]);
00486 }
00487
00488 RETURN(NULL);
00489 }
00490
00491
00492
00493
00494
00495
00496
00497
00498 static int
00499 process_data( maxima_s * M )
00500 {
00501 ENTRY("process_data");
00502 ( void )find_local_maxima( M );
00503
00504 if ( ! create_point_list( M ) )
00505 RETURN(0);
00506
00507 gr_orig_data = M->sdata;
00508 if ( M->negatives )
00509 qsort( M->P.plist, M->P.used, sizeof( int ), point_comp_neg );
00510 else
00511 qsort( M->P.plist, M->P.used, sizeof( int ), point_comp_pos );
00512
00513 if ( M->debug > 1 )
00514 show_point_list_s( "+d point list sorted: ", &M->P, M->debug );
00515
00516 if ( ( M->min_dist > 1.0 ) && ! apply_min_dist( M ) )
00517 RETURN(0);
00518
00519 if ( M->debug > 1 )
00520 show_point_list_s( "+d point list cleaned: ", &M->P, M->debug );
00521
00522 if ( M->outfile )
00523 apply_fill_radius( M );
00524
00525 RETURN(1);
00526 }
00527
00528
00529
00530
00531
00532
00533 static void
00534 show_point_list_s( char * mesg, point_list_s * p, int debug )
00535 {
00536 int c;
00537
00538 ENTRY("show_point_list_s");
00539
00540 if ( mesg ) fputs( mesg, stderr );
00541
00542 fprintf(stderr, "point_list_s @ %p, used = %d, M = %d\n",
00543 (void *)p, p->used, p->M);
00544
00545 if ( debug <= 0 ) EXRETURN;
00546
00547 fprintf(stderr," plist starting @ %p:", (void *)p->plist );
00548
00549 for ( c = 0; c < p->used; c++ )
00550 fprintf(stderr," %d", p->plist[c] );
00551 fprintf(stderr,"\n");
00552
00553 EXRETURN;
00554 }
00555
00556
00557
00558
00559
00560
00561
00562
00563 static void
00564 show_maxima_s( char * mesg, maxima_s * M )
00565 {
00566 ENTRY("show_maxima_s");
00567
00568 if ( mesg ) fputs( mesg, stderr );
00569
00570 fprintf( stderr,
00571 "------------------------------\n"
00572 "dset * : %p\n"
00573 "sdata * : %p\n"
00574 "result * : %p\n"
00575 "nx : %d\n"
00576 "ny : %d\n"
00577 "nz : %d\n"
00578 "nxy : %d\n"
00579 "nvox : %d\n"
00580
00581 "P.plist : %p\n"
00582 "P.used : %d\n"
00583 "P.M : %d\n"
00584
00585 "extrema count : %d\n"
00586
00587 "data_type : %d\n"
00588 "adn_type : %d\n"
00589 "func_type : %d\n"
00590
00591 "outfile : %s\n"
00592 "sval_style : %d\n"
00593
00594 "cutoff : %f\n"
00595 "min_dist : %f\n"
00596 "out_rad : %f\n"
00597
00598 "negatives : %d\n"
00599 "ngbr_style : %d\n"
00600 "overwrite : %d\n"
00601 "quiet : %d\n"
00602 "true_max : %d\n"
00603 "dicom_coords : %d\n"
00604 "debug : %d\n"
00605 "------------------------------\n",
00606
00607 (void *)M->dset, (void *)M->sdata, (void *)M->result,
00608 M->nx, M->ny, M->nz, M->nxy, M->nvox,
00609 (void *)M->P.plist, M->P.used, M->P.M,
00610 M->extrema_count,
00611 M->data_type, M->adn_type, M->func_type,
00612 M->outfile, M->sval_style,
00613 M->cutoff, M->min_dist, M->out_rad,
00614 M->negatives, M->ngbr_style, M->overwrite,
00615 M->quiet, M->true_max, M->dicom_coords, M->debug
00616 );
00617
00618 EXRETURN;
00619 }
00620
00621
00622
00623
00624
00625
00626
00627
00628 static int
00629 apply_min_dist( maxima_s * M )
00630 {
00631 int * iptr, count;
00632 point_list_s newP = { NULL, 0, 0 };
00633
00634
00635 for ( count = 0, iptr = M->P.plist; count < M->P.used; count++, iptr++ )
00636 clear_around_point( *iptr, M, &newP );
00637
00638 free( M->P.plist );
00639 M->P = newP;
00640
00641 return 1;
00642 }
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653 static int
00654 clear_around_point( int p, maxima_s * M, point_list_s * newP )
00655 {
00656 int X, Y, Z;
00657 int xmin, xmax, ymin, ymax, zmin, zmax;
00658 int yc, zc, xrad, yrad, yrad2;
00659 int xbase, ybase, zbase;
00660
00661 short * optr;
00662 float radius = M->min_dist;
00663
00664 static point_list_s P = { NULL, 0, 0 };
00665 P.used = 0;
00666
00667
00668 X = p % M->nx;
00669 Y = (p % M->nxy) / M->nx;
00670 Z = p / M->nxy;
00671
00672 if ( ! *(M->result + ( Z*M->ny + Y ) * M->nx + X ) )
00673 return 1;
00674
00675 zmin = ( Z < radius ) ? Z : radius;
00676 zmax = ( Z + radius >= M->nz ) ? ( M->nz-Z-1 ) : radius;
00677
00678 if ( M->debug > 1 )
00679 fprintf(stderr,"+d index %d, val %f\n", p, M->sdata[p]*gr_fac);
00680
00681 for ( zc = -zmin; zc <= zmax; zc++ )
00682 {
00683 zbase = ( Z + zc ) * M->nx * M->ny;
00684 yrad2 = radius * radius - zc * zc;
00685 yrad = (int)sqrt( yrad2 );
00686
00687 ymin = ( Y < yrad ) ? Y : yrad;
00688 ymax = ( Y + yrad >= M->ny ) ? ( M->ny - Y - 1 ) : yrad;
00689
00690 for ( yc = -ymin; yc <= ymax; yc++ )
00691 {
00692 ybase = ( Y + yc ) * M->nx;
00693 xrad = (int)sqrt( yrad2 - yc * yc );
00694
00695 xmin = ( X < xrad ) ? X : xrad;
00696 xmax = ( X + xrad >= M->nx ) ? ( M->nx - X - 1 ) : xrad;
00697
00698 optr = M->result + ybase + zbase;
00699
00700 for ( xbase = X-xmin; xbase <= X+xmax; xbase++ )
00701 if ( *( optr + xbase ) )
00702 {
00703 *(optr + xbase) = 0;
00704
00705 if ( M->ngbr_style == MAX_WEIGHTED_AVE_STYLE )
00706 {
00707 if ( ! add_point_to_list( &P, xbase+ybase+zbase ) )
00708 return 0;
00709 if ( M->debug > 2 )
00710 fprintf(stderr," coords %d %d %d [%d], value %f\n",
00711 xbase, Y+yc, Z+zc, xbase+ybase+zbase,
00712 M->sdata[xbase+ybase+zbase]*gr_fac);
00713 }
00714 }
00715 }
00716 }
00717
00718 switch( M->ngbr_style )
00719 {
00720 case MAX_SORT_N_REMOVE_STYLE :
00721 if ( ! add_point_to_list( newP, p ) )
00722 return 0;
00723
00724 break;
00725
00726 case MAX_WEIGHTED_AVE_STYLE :
00727 if ( ( p = weighted_index( &P, M ) ) < 0 )
00728 return 0;
00729
00730 if ( ! add_point_to_list( newP, p ) )
00731 return 0;
00732
00733 break;
00734
00735 default :
00736 sprintf( grMessage, "Error: cap_00\nInvalid ngbr_style %d.",
00737 M->ngbr_style );
00738 rERROR( grMessage );
00739 return 0;
00740
00741 break;
00742 }
00743
00744 return 1;
00745 }
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755 static int
00756 weighted_index( point_list_s * P, maxima_s * M )
00757 {
00758 double total_x = 0.0, total_y = 0.0, total_z = 0.0;
00759 double weight = 0.0;
00760 double value;
00761 int x, y, z;
00762 int count, index;
00763 int * iptr;
00764
00765 if ( ( P->plist == NULL ) || ( P->used <= 0 ) )
00766 {
00767 rERROR( "Error wi_00\nEmpty point list." );
00768 return( -1 );
00769 }
00770
00771 if ( P->used == 1 )
00772 return( P->plist[0] );
00773
00774 for ( count = 0, iptr = P->plist; count < P->used; count++, iptr++ )
00775 {
00776 index = *iptr;
00777
00778 x = index % M->nx;
00779 y = ( index % M->nxy ) / M->nx;
00780 z = index / M->nxy;
00781
00782 value = M->sdata[ index ];
00783
00784 weight += value;
00785 total_x += value * x;
00786 total_y += value * y;
00787 total_z += value * z;
00788 }
00789
00790 if ( M->debug > 1 )
00791 fprintf(stderr, "-d nvals, weight, ave = %d, %f, %f\n",
00792 P->used, weight*gr_fac, weight*gr_fac/P->used);
00793
00794 if ( weight <= 0.0 )
00795 {
00796 sprintf( grMessage, "Error: wi_10\nunexpected weight of %f", weight );
00797 rERROR( grMessage );
00798 }
00799
00800 x = ( int )( total_x / weight + 0.4 );
00801 y = ( int )( total_y / weight + 0.4 );
00802 z = ( int )( total_z / weight + 0.4 );
00803
00804 index = ( z * M->ny + y ) * M->nx + x;
00805
00806 if ( M->debug > 1 )
00807 fprintf(stderr, "-d weighted i,j,k, ind, val = %d, %d, %d, %d, %f\n",
00808 x, y, z, index, M->sdata[index]*gr_fac);
00809
00810 return index;
00811 }
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822 static int
00823 create_point_list( maxima_s * M )
00824 {
00825 short * mptr;
00826 int count;
00827 point_list_s * P = &M->P;
00828
00829 ENTRY("create_pint_list");
00830
00831 mptr = M->result;
00832 for ( count = 0; count < M->nvox; count++ )
00833 if ( *mptr++ )
00834 if ( ! add_point_to_list( P, count ) )
00835 RETURN(0);
00836
00837 if ( M->debug > 0 )
00838 show_point_list_s( "+d point list created: ", &M->P, M->debug );
00839
00840 RETURN(1);
00841 }
00842
00843
00844
00845
00846
00847
00848
00849
00850 static int
00851 add_point_to_list( point_list_s * P, int offset )
00852 {
00853 ENTRY("add_point_to_list");
00854 if ( ! P->plist )
00855 {
00856 P->M = 100;
00857 P->used = 0;
00858
00859 if ( ( P->plist = (int *)malloc( P->M * sizeof(int) ) ) == NULL )
00860 {
00861 rERROR( "Error: aptl_10\n"
00862 "Failed to allocate memory for initial point list.\n" );
00863 RETURN(0);
00864 }
00865 }
00866 else if ( P->used == P->M )
00867 {
00868 P->M = P->M + 100;
00869
00870 if ( ( P->plist = (int *)realloc( P->plist, P->M*sizeof(int))) == NULL )
00871 {
00872 sprintf( grMessage, "Error: aptl_20\n"
00873 "Failed to reallocate %d ints for point list", P->M );
00874 RETURN(0);
00875 }
00876 }
00877
00878 P->plist[P->used] = offset;
00879 P->used++;
00880
00881 RETURN(1);
00882 }
00883
00884
00885
00886
00887
00888
00889
00890
00891 static int
00892 apply_fill_radius( maxima_s * M )
00893 {
00894 int count, outval;
00895 int x, y, z;
00896 int * iptr;
00897
00898 ENTRY("apply_fill_radius");
00899
00900 for ( count = 0, iptr = M->P.plist; count < M->P.used; count++, iptr++ )
00901 {
00902 outval = (M->sval_style == 0) ? MAX_MASK_FILL_VAL :
00903 (M->sval_style == 1) ? (count+1) :
00904 (M->P.used - count);
00905
00906 if ( M->out_rad < 1.0 )
00907 {
00908 M->result[ *iptr ] = outval;
00909 continue;
00910 }
00911
00912 x = *iptr % M->nx;
00913 y = ( *iptr % M->nxy ) / M->nx;
00914 z = *iptr / M->nxy;
00915
00916 radial_fill( x, y, z, M, outval );
00917 }
00918
00919 RETURN(1);
00920 }
00921
00922
00923
00924
00925
00926
00927
00928
00929 static int
00930 radial_fill( int X, int Y, int Z, maxima_s * M, int val )
00931 {
00932 int xmin, xmax, ymin, ymax, zmin, zmax;
00933 int yc, zc, xrad, yrad, yrad2;
00934 int xbase, ybase, zbase;
00935
00936 short * sptr, * optr;
00937 float radius = M->out_rad;
00938
00939 ENTRY("radial_fill");
00940
00941 zmin = ( Z < radius ) ? Z : radius;
00942 zmax = ( Z + radius >= M->nz ) ? ( M->nz-Z-1 ) : radius;
00943
00944 for ( zc = -zmin; zc <= zmax; zc++ )
00945 {
00946 zbase = ( Z + zc ) * M->nx * M->ny;
00947 yrad2 = radius * radius - zc * zc;
00948 yrad = (int)sqrt( yrad2 );
00949
00950 ymin = ( Y < yrad ) ? Y : yrad;
00951 ymax = ( Y + yrad >= M->ny ) ? ( M->ny - Y - 1 ) : yrad;
00952
00953 for ( yc = -ymin; yc <= ymax; yc++ )
00954 {
00955 ybase = ( Y + yc ) * M->nx;
00956 xrad = (int)sqrt( yrad2 - yc * yc );
00957
00958 xmin = ( X < xrad ) ? X : xrad;
00959 xmax = ( X + xrad >= M->nx ) ? ( M->nx - X - 1 ) : xrad;
00960
00961 optr = M->result + ybase + zbase;
00962
00963 for ( xbase = X-xmin; xbase <= X+xmax; xbase++ )
00964 {
00965 sptr = optr + xbase;
00966
00967 if ( ! *sptr )
00968 *sptr = val;
00969 }
00970 }
00971 }
00972
00973 RETURN(1);
00974 }
00975
00976
00977
00978
00979
00980
00981
00982
00983 static int
00984 find_local_maxima( maxima_s * M )
00985 {
00986 short * sourcep, * destp;
00987 int c, cx, cy, cz;
00988 int maxx = M->nx - 1;
00989 int maxy = M->ny - 1;
00990 int nx = M->nx, nxy = M->nx * M->ny;
00991 int offset[ 26 ];
00992
00993 ENTRY("find_local_maxima");
00994
00995 offset[ 0] = +1;
00996 offset[ 1] = -1;
00997 offset[ 2] = nx;
00998 offset[ 3] = nx+1;
00999 offset[ 4] = nx-1;
01000 offset[ 5] = -nx;
01001 offset[ 6] = -nx+1;
01002 offset[ 7] = -nx-1;
01003 offset[ 8] = nxy;
01004 offset[ 9] = nxy+1;
01005 offset[10] = nxy-1;
01006 offset[11] = nxy+nx;
01007 offset[12] = nxy+nx+1;
01008 offset[13] = nxy+nx-1;
01009 offset[14] = nxy-nx;
01010 offset[15] = nxy-nx+1;
01011 offset[16] = nxy-nx-1;
01012 offset[17] = -nxy;
01013 offset[18] = -nxy+1;
01014 offset[19] = -nxy-1;
01015 offset[20] = -nxy+nx;
01016 offset[21] = -nxy+nx+1;
01017 offset[22] = -nxy+nx-1;
01018 offset[23] = -nxy-nx;
01019 offset[24] = -nxy-nx+1;
01020 offset[25] = -nxy-nx-1;
01021
01022 sourcep = M->sdata + nxy;
01023 destp = M->result + nxy;
01024
01025 for ( cz = 0; cz < M->nz-2; cz++ )
01026 {
01027 for ( cy = 0; cy < M->ny; cy++ )
01028 {
01029 if ( ( cy == 0 ) || ( cy == maxy ) )
01030 {
01031 sourcep += nx;
01032 destp += nx;
01033
01034 continue;
01035 }
01036
01037 for ( cx = 0; cx < M->nx; cx++ )
01038 {
01039 if ( ( cx == 0 ) || ( cx == maxx ) )
01040 {
01041 sourcep++;
01042 destp++;
01043
01044 continue;
01045 }
01046
01047 if ( ! M->negatives && ( *sourcep < M->cutoff ) )
01048 {
01049 sourcep++;
01050 destp++;
01051
01052 continue;
01053 }
01054
01055 if ( M->negatives && ( *sourcep > M->cutoff ) )
01056 {
01057 sourcep++;
01058 destp++;
01059
01060 continue;
01061 }
01062
01063 *destp = MAX_MASK_FILL_VAL;
01064 M->extrema_count++;
01065
01066 if ( M->true_max )
01067 {
01068 if ( ! M->negatives )
01069 {
01070 for ( c = 0; c < 26; c++ )
01071 if ( *sourcep <= sourcep[offset[c]] )
01072 {
01073 *destp = 0;
01074 M->extrema_count--;
01075
01076 break;
01077 }
01078 }
01079 else
01080 {
01081 for ( c = 0; c < 26; c++ )
01082 if ( *sourcep >= sourcep[offset[c]] )
01083 {
01084 *destp = 0;
01085 M->extrema_count--;
01086
01087 break;
01088 }
01089 }
01090 }
01091 else
01092 {
01093 if ( ! M->negatives )
01094 {
01095 for ( c = 0; c < 26; c++ )
01096 if ( *sourcep < sourcep[offset[c]] )
01097 {
01098 *destp = 0;
01099 M->extrema_count--;
01100
01101 break;
01102 }
01103 }
01104 else
01105 {
01106 for ( c = 0; c < 26; c++ )
01107 if ( *sourcep > sourcep[offset[c]] )
01108 {
01109 *destp = 0;
01110 M->extrema_count--;
01111
01112 break;
01113 }
01114 }
01115 }
01116
01117 sourcep++;
01118 destp++;
01119 }
01120 }
01121 }
01122
01123 if ( M->debug > 3 ) show_maxima_s( "post find local maxima: ", M );
01124
01125 RETURN(1);
01126 }
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137 static int
01138 write_results( r_afni_s * A, maxima_s * M, PLUGIN_interface * plint )
01139 {
01140 THD_3dim_dataset * newdset;
01141
01142 ENTRY("write_results");
01143
01144 if ( ! M->quiet )
01145 display_coords( A, M );
01146
01147 if ( ! *M->outfile )
01148 RETURN(1);
01149
01150
01151
01152
01153 if ( ( newdset = EDIT_empty_copy( M->dset ) ) == NULL )
01154 {
01155 rERROR( "Error: wr_00\n" "Failed to copy dataset." );
01156 RETURN(0);
01157 }
01158
01159 { char * his = PLUTO_commandstring(plint) ;
01160 tross_Copy_History( M->dset , newdset ) ;
01161 tross_Append_History( newdset , his ) ; free(his) ;
01162 }
01163
01164 EDIT_dset_items( newdset,
01165 ADN_prefix, M->outfile,
01166 ADN_label1, M->outfile,
01167 ADN_nvals, 1,
01168 ADN_ntt, 0,
01169 ADN_type, HEAD_FUNC_TYPE,
01170 ADN_func_type, FUNC_FIM_TYPE,
01171 ADN_none
01172 );
01173
01174 EDIT_substitute_brick( newdset, 0, M->data_type, M->result );
01175 EDIT_BRICK_FACTOR ( newdset, 0, 0.0 );
01176
01177 if ( PLUTO_add_dset( plint, newdset, DSET_ACTION_MAKE_CURRENT ) )
01178 {
01179 rERROR( "Error: wr_10\n" "Failed to make current dataset." );
01180 RETURN(0);
01181 }
01182 else
01183 DSET_unload( M->dset );
01184
01185 RETURN(1);
01186 }
01187
01188
01189
01190
01191
01192
01193
01194
01195 static int
01196 display_coords( r_afni_s * A, maxima_s * M )
01197 {
01198 THD_fvec3 f3;
01199 THD_ivec3 i3;
01200 float prod, factor = A->factor[0];
01201 short * optr;
01202 short * mptr;
01203 int * iptr;
01204 int X, Y, Z, count;
01205
01206 point_list_s * P = &M->P;
01207
01208 ENTRY("display_coords");
01209
01210 printf( "---------------------------------------------\n" );
01211 if ( M->dicom_coords )
01212 printf( "RAI mm coordinates:\n\n" );
01213 else
01214 printf( "%c%c%c mm coordinates:\n\n",
01215 ORIENT_typestr[M->dset->daxes->xxorient][0],
01216 ORIENT_typestr[M->dset->daxes->yyorient][0],
01217 ORIENT_typestr[M->dset->daxes->zzorient][0] );
01218
01219 for ( count = 0, iptr = P->plist; count < P->used; count++, iptr++ )
01220 {
01221 X = *iptr % M->nx;
01222 Y = (*iptr % M->nxy) / M->nx;
01223 Z = *iptr / M->nxy;
01224 i3.ijk[0] = X; i3.ijk[1] = Y; i3.ijk[2] = Z;
01225 f3 = THD_3dind_to_3dmm_no_wod(M->dset, i3);
01226 if ( M->dicom_coords )
01227 f3 = THD_3dmm_to_dicomm(M->dset, f3);
01228
01229 optr = M->sdata + *iptr;
01230 mptr = M->result + *iptr;
01231
01232 if ( factor == 1 )
01233 {
01234
01235 printf( "(%7.2f %7.2f %7.2f) : val = %d\n",
01236 f3.xyz[0], f3.xyz[1], f3.xyz[2], *optr );
01237 }
01238 else
01239 {
01240 prod = *optr * factor;
01241
01242 printf( "(%7.2f %7.2f %7.2f) : val = %f\n",
01243 f3.xyz[0], f3.xyz[1], f3.xyz[2], prod );
01244 }
01245 }
01246
01247 if ( P->used )
01248 printf( "\nnumber of extrema = %d\n", P->used );
01249 else
01250 printf( "No extrema found.\n" );
01251 printf( "---------------------------------------------\n" );
01252
01253
01254 RETURN(1);
01255 }
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266 static int
01267 init_afni_s( r_afni_s * A )
01268 {
01269 ENTRY("init_afni_s");
01270
01271 memset( A, 0, sizeof( r_afni_s ) );
01272
01273 A->must_be_short = 1;
01274 A->want_floats = 1;
01275 A->subs_must_equal = 1;
01276 A->max_subs = 1;
01277
01278 RETURN(1);
01279 }
01280
01281
01282
01283
01284
01285
01286
01287
01288 static int
01289 init_maxima_s( maxima_s * M, r_afni_s * A, char * outprefix )
01290 {
01291 ENTRY("init_maxima_s");
01292
01293 M->dset = A->dset[0];
01294
01295 M->sdata = A->simage[0];
01296
01297 if ( ( M->result = (short *)calloc( A->nvox, sizeof( short ) ) ) == NULL )
01298 {
01299 sprintf( grMessage, "Error: ims_05\n"
01300 "Failed to allocate M for %d shorts.", A->nvox );
01301 rERROR( grMessage );
01302 RETURN(0);
01303 }
01304
01305 M->nx = A->nx;
01306 M->ny = A->ny;
01307 M->nz = A->nz;
01308 M->nxy = A->nx * A->ny;
01309 M->nvox = A->nvox;
01310
01311 M->P.plist = NULL;
01312 M->P.used = 0;
01313 M->P.M = 0;
01314
01315 M->extrema_count = 0;
01316
01317 M->data_type = MRI_short;
01318 M->adn_type = HEAD_FUNC_TYPE;
01319 M->func_type = FUNC_FIM_TYPE;
01320
01321 if ( outprefix && strlen( outprefix ) > R_FILE_L )
01322 {
01323 sprintf( grMessage, "Error: ims_10\n"
01324 "Outfile prefix exceeds %d characters.", R_FILE_L );
01325 rERROR( grMessage );
01326 RETURN(0);
01327 }
01328
01329 if ( outprefix )
01330 strcpy( M->outfile, outprefix );
01331 else
01332 *M->outfile = 0;
01333
01334 M->cutoff = 0.0;
01335 M->min_dist = 0.0;
01336 M->out_rad = 0.0;
01337
01338 M->negatives = 0;
01339 M->ngbr_style = MAX_SORT_N_REMOVE_STYLE;
01340 M->overwrite = 0;
01341 M->quiet = 0;
01342 M->true_max = 0;
01343 M->debug = 0;
01344
01345 RETURN(1);
01346 }
01347
01348
01349
01350
01351
01352
01353
01354
01355 int
01356 point_comp_pos( const void * p1, const void * p2 )
01357 {
01358 short v1 = *( gr_orig_data + *(int *)p1 );
01359 short v2 = *( gr_orig_data + *(int *)p2 );
01360
01361 if ( v1 < v2 )
01362 return 1;
01363 else if ( v1 > v2 )
01364 return -1;
01365
01366 return 0;
01367 }
01368
01369
01370
01371
01372
01373
01374
01375
01376 int
01377 point_comp_neg( const void * p1, const void * p2 )
01378 {
01379 short v1 = *( gr_orig_data + *(int *)p1 );
01380 short v2 = *( gr_orig_data + *(int *)p2 );
01381
01382
01383 if ( v1 > v2 )
01384 return 1;
01385 else if ( v1 < v2 )
01386 return -1;
01387
01388
01389 return 0;
01390 }
01391
01392
01393
01394
01395
01396
01397
01398
01399 int
01400 r_set_afni_s_from_dset( r_afni_s * A, THD_3dim_dataset * dset )
01401 {
01402 ENTRY("r_set_afni_s_from_dset");
01403
01404 if ( A->num_dsets >= R_MAX_AFNI_DSETS )
01405 {
01406 sprintf( grMessage, "Error: rsasfd_00\n"
01407 "We only have memory to hold %d datasets. exiting...\n",
01408 R_MAX_AFNI_DSETS );
01409 rERROR( grMessage );
01410
01411 RETURN(0);
01412 }
01413
01414 A->dset[ 0 ] = dset;
01415 A->simage[ 0 ] = ( short * )DSET_ARRAY( dset, A->sub_brick );
01416
01417 if ( !A->simage[0] )
01418 {
01419 sprintf(grMessage,
01420 "** data not available, is this in warp-on-demand mode?\n");
01421 rERROR(grMessage);
01422 RETURN(0);
01423 }
01424
01425 if ((A->factor[0] = DSET_BRICK_FACTOR(dset, A->sub_brick)) == 0.0 )
01426 A->factor[0] = 1.0;
01427
01428 A->subs [ 0 ] = DSET_NVALS( dset );
01429
01430 A->nx = dset->daxes->nxx;
01431 A->ny = dset->daxes->nyy;
01432 A->nz = dset->daxes->nzz;
01433 A->nvox = A->nx * A->ny * A->nz;
01434
01435 if ( A->want_floats )
01436 {
01437 int count;
01438 short * sptr;
01439 float * fptr;
01440 float factor = A->factor[ 0 ];
01441
01442 if ( ( A->fimage[ 0 ] =
01443 ( float * )calloc( A->nvox, sizeof( float ) ) ) == NULL )
01444 {
01445 sprintf( grMessage, "Error: rsasfd_10\n"
01446 "Failed to allocate memory for %d floats.\n", A->nvox );
01447 rERROR( grMessage );
01448
01449 RETURN(0);
01450 }
01451
01452 fptr = A->fimage[ 0 ];
01453 sptr = A->simage[ 0 ];
01454 for ( count = 0; count < A->nvox; count++ )
01455 *fptr++ = *sptr++ * factor;
01456 }
01457
01458 A->max_u_short = r_get_max_u_short( (u_short *)A->simage[0], A->nvox );
01459
01460
01461
01462 RETURN(1);
01463 }
01464
01465
01466
01467
01468
01469
01470
01471
01472 u_short
01473 r_get_max_u_short( u_short * S, int size )
01474 {
01475 u_short * usptr, max = *S;
01476 int c = 0;
01477
01478 for ( c = 0, usptr = S; c < size; c++, usptr++ )
01479 {
01480 if ( *usptr > max )
01481 max = *usptr;
01482 }
01483
01484 return max;
01485 }
01486
01487
01488
01489
01490
01491
01492
01493
01494 static void
01495 free_memory( r_afni_s * A, maxima_s * M )
01496 {
01497 ENTRY("free_memory");
01498 if ( A->want_floats && A->fimage[0] )
01499 free( A->fimage[0] );
01500
01501 if ( M->result && !M->outfile[0] )
01502 free( M->result );
01503
01504 if ( M->P.plist )
01505 free( M->P.plist );
01506
01507 EXRETURN;
01508 }
01509
01510
