1 /*****************************************************************************
2    Major portions of this software are copyrighted by the Medical College
3    of Wisconsin, 1994-2000, and are released under the Gnu General Public
4    License, Version 2.  See the file README.Copyright for details.
5 ******************************************************************************/
6 
7 /*----------------------------------------------------------------------
8  *
9  *  plug_maxima.c	- AFNI plugin to locate relative extrema
10  *
11  *----------------------------------------------------------------------
12 */
13 
14 /*----------------------------------------------------------------------
15  * history:  the history is maintained in 3dmaxima
16  *----------------------------------------------------------------------
17 */
18 
19 #include "afni.h"
20 #include "maxima.h"
21 
22 #ifndef ALLOW_PLUGINS
23 #  error "Plugins not properly set up -- see machdep.h"
24 #endif
25 
26 /***********************************************************************
27   plugin to find local extrema
28 ************************************************************************/
29 
30 char * MAXIMA_main( PLUGIN_interface * );
31 char * process_args( r_afni_s * A, maxima_s * M, PLUGIN_interface * plint );
32 
33 static char      * grStyle[] = { "Sort-n-Remove", "Weighted-Average" };
34 static char      * grSvals[] = { "1 (default)", "1 to N", "N to 1" };
35 static char      * grNY[]    = { "No", "Yes" };
36 
37 static char        helpstring[] =
38     "Maxima - used to locate extrema in a functional dataset.\n"
39     "\n"
40     "This plugin reads a functional dataset and locates any relative extrema\n"
41     "(maximums or minimums, depending on the user option).  A _relative_\n"
42     "maximum is a point that is greater than all neighbors (not necessarily\n"
43     "greater than all other values in the sub-brick).  The output from this\n"
44     "process can be text based (sent to the terminal window) and it can be a\n"
45     "mask (integral) dataset, where the locations of the extrema are set.\n"
46     "\n"
47     "When writing a dataset, it is often useful to set a sphere around each\n"
48     "extrema, not to just set individual voxels.  This makes viewing those\n"
49     "locations much more reasonable.  Also, if the 'Sphere Values' option is\n"
50     "set to 'N to 1', the sphere around the most extreme voxel will get the\n"
51     "value N, giving it the 'top' color in afni (and so on, down to 1).\n"
52     "\n"
53     "Notes : The only required option is the input dataset.\n"
54     "        Input datasets must be of type short.\n"
55     "        All distances are in voxel units.\n"
56     "\n"
57     "                        ***  Options  ***\n"
58     "\n"
59     "-----  Input:  -----\n"
60     "\n"
61     "   Dataset   - dataset to locate extrema within\n"
62     "\n"
63     "   Sub-brick - sub-brick index of volume to locate extrema within\n"
64     "\n"
65     "-----  Output Dset:  -----\n"
66     "\n"
67     "   Prefix    - prefix for an output functional dataset\n"
68     "\n"
69     "               This dataset may be viewed as a mask.  Its set points\n"
70     "               will correspond to any selected extrema.  The \"Output\n"
71     "               Size\" option will change those points to spheres.\n"
72     "\n"
73     "   Values    - values to fill output spheres with\n"
74     "\n"
75     "               1 (default) - fill all spheres with the value 1\n"
76     "               1 to N      - fill sphere of most extreme point with 1,\n"
77     "                             the next with 2, and so on\n"
78     "               N to 1      - fill sphere of most extreme point with N,\n"
79     "                             the next with N-1, etc., perhaps matching\n"
80     "                             Overlay coloring better\n"
81     "\n"
82     "-----  Threshold:  -----\n"
83     "\n"
84     "   Cutoff    - provides a cutoff value for extrema\n"
85     "\n"
86     "               Extrema not meeting this cutoff will not be considered.\n"
87     "               Note that if the \"Neg Extrema\" option is applied, the\n"
88     "               user will generally want a negative threshold.\n"
89     "\n"
90     "-----  Separation:  -----\n"
91     "\n"
92     "   Distance  - minimum acceptable distance between extrema\n"
93     "\n"
94     "               Less significant extrema which are too close to more\n"
95     "               significant extrema will be discounted in some way.\n"
96     "               This option works with the \"Neighbor Style\" option.\n"
97     "\n"
98     "               Note that the distance is in voxels, not mm.\n"
99     "\n"
100     "-----  Output Size:  -----\n"
101     "\n"
102     "   Radius    - output mask radius around set points\n"
103     "\n"
104     "               If the user wants the output BRIK to consist of spheres\n"
105     "               centered around extrema points, this option can be used\n"
106     "               to set the radius for those spheres.\n"
107     "\n"
108     "-----  Neighbor:  -----\n"
109     "\n"
110     "   Style     - technique for handling clustered extrema\n"
111     "\n"
112     "               If extrema are not as far apart as is specified by the\n"
113     "               \"Separation\" option, the Neighbor (Style) option\n"
114     "               specifies how to handle those points.\n"
115     "\n"
116     "               Sort-n-Remove : The extrema are sorted by magnitude.\n"
117     "                  For each such extrema, if the extrema is still set\n"
118     "                  (not previously removed), all (less significant)\n"
119     "                  neighbors within the Separation radius are removed.\n"
120     "\n"
121     "               Weighted-Average : Again, traverse the sorted list of\n"
122     "                  extrema.  Set the current extrema to the center of\n"
123     "                  mass of all extrema within the Separation radius.\n"
124     "\n"
125     "-----  Params:  -----\n"
126     "\n"
127     "   Neg Extr  - search for negative extrema (minima)\n"
128     "\n"
129     "               This will search for the minima of the dataset.\n"
130     "               Note that a negative threshold may be desired.\n"
131     "\n"
132     "   True Max  - do not consider extrema with an equal valued neighbor\n"
133     "\n"
134     "               As a default points may be considered extrema, even if\n"
135     "               they have a neighbor of the same value.  Setting this\n"
136     "               option requires extrema to be strictly greater than\n"
137     "               any of their neighbors.\n"
138     "               Note that with this set, extrema locations that have\n"
139     "               more than one voxel at the maximal value are ignored.\n"
140     "\n"
141     "-----  Output Text:  -----\n"
142     "\n"
143     "   No Text  - do not display the extrma points as text\n"
144     "\n"
145     "   Debug    - how much extra information to output to the terminal\n"
146     "\n"
147     "-----  Output Coords:  -----\n"
148     "\n"
149     "   Dicom    - Yes/No: whether to display output in Dicom orientation\n"
150     "\n"
151     "              In Dicom orientation, the output is RAI, meaning right,\n"
152     "              anterior and inferior are negative coordinates, and they\n"
153     "              are printed in that order (RL, then AP, then IS).\n"
154     "\n"
155     "              If this option is set to NO, the dataset orientation is\n"
156     "              used, whichever of the 48 it happens to be.\n"
157     "\n"
158     "              Note that in either case, the output orientation is\n"
159     "              printed above the results, in the terminal window.\n"
160     "\n"
161     "Author: (the hopefully-no-longer-lamented) R Reynolds\n"
162     "\n"
163    ;
164 
165 
166 
167 /***********************************************************************
168    Set up the interface to the user
169 ************************************************************************/
170 
171 
172 DEFINE_PLUGIN_PROTOTYPE
173 
PLUGIN_init(int ncall)174 PLUGIN_interface * PLUGIN_init( int ncall )
175 {
176    PLUGIN_interface * plint ;
177 
178    if( ncall > 0 ) return NULL ;  /* only one interface */
179 
180    /* create the new interface */
181 
182    plint = PLUTO_new_interface( "Maxima", "find extrema in a dataset",
183 		helpstring, PLUGIN_CALL_VIA_MENU , MAXIMA_main );
184 
185    PLUTO_add_hint( plint, "find local maxima/minima" );
186 
187    PLUTO_set_sequence( plint , "z:Reynolds" ) ;
188 
189    /*-- first line of input: input dataset --*/
190 
191    PLUTO_add_option( plint, "Input" , "Input" , TRUE );
192    PLUTO_add_hint( plint, "choose dataset for input" );
193    PLUTO_add_dataset(plint, "Dataset" , ANAT_ALL_MASK , FUNC_ALL_MASK,
194                                          DIMEN_ALL_MASK | BRICK_SHORT_MASK );
195    PLUTO_add_number( plint , "Sub-brick" , 0,9999,0 , 0,1 ) ; /* new [rickr] */
196 
197 
198    /*-- second line of input: prefix for output dataset --*/
199 
200    PLUTO_add_option( plint, "Output Dset" , "prefix" , FALSE );
201    PLUTO_add_hint( plint, "options for the creation of an output dataset");
202    PLUTO_add_string( plint, "Prefix", 0 , NULL, 19 );
203    PLUTO_add_hint( plint, "option: choose dataset prefix for output" );
204    PLUTO_add_string( plint, "Sphere Values", 3 , grSvals, 0 );
205    PLUTO_add_hint( plint, "option: choose value style for output spheres" );
206 
207    /*-- third line of input: cutoff option --*/
208 
209    PLUTO_add_option( plint, "Threshold" , "cutoff" , FALSE ) ;
210    PLUTO_add_hint( plint, "option: choose a threshold for value at extrema" );
211    PLUTO_add_number( plint, "Cutoff", 0, 0, 0, 1000, 1 );
212 
213    /*-- fourth line of input: min_dist option --*/
214 
215    PLUTO_add_option( plint, "Separation" , "min_dist" , FALSE ) ;
216    PLUTO_add_hint( plint, "option: choose a minimum distance between extrema" );
217    PLUTO_add_number( plint, "Distance(vox)", 0, 1000, 1, 40, 1 );
218 
219    /*-- fifth line of input: out_rad option --*/
220 
221    PLUTO_add_option( plint, "Output Size" , "out_rad" , FALSE ) ;
222    PLUTO_add_hint( plint, "option: choose a spherical radius around extrema "
223 			  "points in mask" );
224    PLUTO_add_number( plint, "Radius(vox)", 0, 1000, 1, 50, 1 );
225 
226    /*-- sixth line of input: style option --*/
227 
228    PLUTO_add_option( plint, "Neighbor" , "style" , FALSE ) ;
229    PLUTO_add_hint( plint, "option: technique for neighbor removal" );
230    PLUTO_add_string( plint, "Style", 2, grStyle, 0 );
231 
232    /*-- seventh line of input: negatives and true max options --*/
233 
234    PLUTO_add_option( plint, "Params" , "params" , FALSE ) ;
235    PLUTO_add_hint( plint, "options: negative extrema and true max" );
236    PLUTO_add_string( plint, "Neg Extrema", 2, grNY, 0 );
237    PLUTO_add_hint( plint, "search for negative extrema, not positive" );
238    PLUTO_add_string( plint, "True Max", 2, grNY, 0 );
239    PLUTO_add_hint( plint, "exclude extrema with equal neighbors" );
240 
241    /*-- eighth line of input: true_max option --*/
242 
243    PLUTO_add_option( plint, "Output Text" , "output" , FALSE ) ;
244    PLUTO_add_hint( plint, "options: no output text, debug level" );
245    PLUTO_add_string( plint, "No Text Out", 2, grNY, 0 );
246    PLUTO_add_hint( plint, "do not output extrema as text (to terminal)" );
247    PLUTO_add_number( plint, "Debug Level", 0, 4, 0, 0, 0 );
248    PLUTO_add_hint( plint, "search for negative extrema, not positive" );
249 
250    /*-- ninth line of input: dicom_coords option --*/
251 
252    PLUTO_add_option( plint, "Output Coords" , "dicom_coords" , FALSE ) ;
253    PLUTO_add_hint( plint, "option: output coordinates in Dicom format" );
254    PLUTO_add_string( plint, "Dicom Coords", 2, grNY, 1 );
255 
256    return plint;
257 }
258 
259 
260 /*----------------------------------------------------------------------
261 **
262 **  Main routine for this plugin (will be called from AFNI).
263 **
264 **----------------------------------------------------------------------
265 */
MAXIMA_main(PLUGIN_interface * plint)266 char * MAXIMA_main( PLUGIN_interface * plint )
267 {
268     r_afni_s           A;
269     maxima_s           M;
270     char             * ret_string = NULL;
271     THD_3dim_dataset * dset;
272 
273 
274     if ( ( ret_string = process_args( &A, &M, plint ) ) != NULL )
275 	return ret_string;
276 
277     if ( ! process_data( &M ) )
278 	return  "************************************\n"
279 	        "MAXIMA_main: data processing failure\n"
280 		"************************************";
281 
282     if ( ! M.quiet ) display_coords( &A, &M );
283 
284     if ( *M.outfile )
285     {
286         if ( ( dset = EDIT_empty_copy( M.dset ) ) == NULL )
287         {
288             rERROR( "Error: wr_00\n" "Failed to copy dataset." );
289             return(0);
290         }
291 
292         { /* add dataset history */
293           char * his = PLUTO_commandstring(plint) ;
294           tross_Copy_History( M.dset , dset ) ;
295           tross_Append_History( dset , his ) ; free(his) ;
296         }
297 
298         if ( ! set_results( &A, &M, dset ) )
299             return  "***********************************\n"
300                     "MAXIMA_main: result writing failure\n"
301                     "***********************************";
302 
303         if ( PLUTO_add_dset( plint, dset, DSET_ACTION_MAKE_CURRENT ) )
304         {
305             rERROR( "Error: wr_10\n" "Failed to make current dataset." );
306             return(0);
307         }
308         else
309             DSET_unload( M.dset );
310     }
311 
312 
313     free_memory( &A, &M );
314 
315     return NULL;
316 }
317 
318 
319 /*----------------------------------------------------------------------
320 **
321 **  Process arguments.
322 **
323 **  Return a description string on failure.
324 **
325 **----------------------------------------------------------------------
326 */
process_args(r_afni_s * A,maxima_s * M,PLUGIN_interface * plint)327 char * process_args( r_afni_s * A, maxima_s * M, PLUGIN_interface * plint )
328 {
329     THD_3dim_dataset * dset;
330     MCW_idcode       * idc ;
331     char             * optag, * outfile = NULL, * str;
332     float              cutoff = 0.0, min_dist = 0.0, out_rad = 0.0;
333     int                negatives = 0, quiet = 0, true_max = 0, opcnt = 0;
334     int                val, debug = 0, style = MAX_SORT_N_REMOVE_STYLE, sb;
335     int                sval_style = 0, dicom_coords = 1;
336 
337 ENTRY("process_args");
338     /* get AFNI inputs */
339 
340     if( plint == NULL )
341 	RETURN("----------------------\n"
342                "arguments : NULL input\n"
343                "----------------------");
344 
345     if ( ! init_afni_s( A ) )
346 	RETURN( "------------------------\n"
347 		"arguments : init failure\n"
348 		"------------------------");
349 
350     PLUTO_next_option( plint );
351     idc  = PLUTO_get_idcode( plint );
352     dset = PLUTO_find_dset( idc );
353 
354     if( dset == NULL )
355 	RETURN("-----------------------------\n"
356                "arguments : bad input dataset\n"
357                "-----------------------------");
358 
359     sb = (int)PLUTO_get_number( plint );	/* 2004 Feb 20 [rickr] */
360     if ( sb >= DSET_NVALS(dset) || sb < 0 )
361 	RETURN("--------------------------\n"
362                "arguments : bad sub-brick \n"
363                "--------------------------");
364     A->sub_brick = sb;
365 
366     DSET_load( dset );
367 
368     for ( optag  = PLUTO_get_optiontag( plint );
369 	  optag != NULL;
370 	  optag  = PLUTO_get_optiontag( plint )
371 	)
372     {
373 	if ( ! strcmp( optag, "prefix" ) )
374 	{
375 	    outfile = PLUTO_get_string( plint );
376 	    if ( ! PLUTO_prefix_ok( outfile ) )
377 		RETURN( "-------------------------\n"
378 		        "options : bad file prefix\n"
379 			"-------------------------");
380             sval_style = PLUTO_string_index(PLUTO_get_string(plint),3,grSvals);
381 	}
382 	else if ( ! strcmp( optag, "cutoff" ) )
383 	{
384 	    cutoff = PLUTO_get_number( plint );
385 	}
386 	else if ( ! strcmp( optag, "min_dist" ) )
387 	{
388 	    if ( ( min_dist = PLUTO_get_number( plint ) ) < 0 )
389 		RETURN( "-----------------------------------------\n"
390 			"options : Separation must be non-negative\n"
391 			"-----------------------------------------");
392 	}
393 	else if ( ! strcmp( optag, "out_rad" ) )
394 	{
395 	    if ( ( out_rad = PLUTO_get_number( plint ) ) < 0 )
396 		RETURN( "--------------------------------------------\n"
397 			"options : Output radius must be non-negative\n"
398 			"--------------------------------------------");
399 	}
400 	else if ( ! strcmp( optag, "params" ) )
401 	{
402 	    str = PLUTO_get_string( plint );		/* Neg Extrema */
403 	    val = PLUTO_string_index(str, 2, grNY);
404 	    if ( val > 0 ) negatives = 1;
405 	    str = PLUTO_get_string( plint );		/* True Max    */
406 	    val = PLUTO_string_index(str, 2, grNY);
407 	    if ( val > 0 ) true_max = 1;
408 	}
409 	else if ( ! strcmp( optag, "style" ) )
410 	{
411 	    if ( ( str = PLUTO_get_string( plint ) ) == NULL )
412 		RETURN( "-------------------------------\n"
413 			"options : missing style string?\n"
414 			"-------------------------------");
415 	    if ((( style = PLUTO_string_index(str, 2, grStyle))
416 				< 0 ) || ( style >= MAX_MAX_STYLE ) )
417 	    {
418 		sprintf( grMessage,
419 		    "---------------------------\n"
420 		    "options : bad style is %d\n"
421 		    "---------------------------", style );
422 		RETURN(grMessage);
423 	    }
424 	    style++;
425 	}
426 	else if ( ! strcmp( optag, "output" ) )
427 	{
428 	    str = PLUTO_get_string( plint );		/* No Text Out */
429 	    val = PLUTO_string_index(str, 2, grNY);
430 	    if ( val > 0 ) quiet = 1;
431 	    debug = PLUTO_get_number( plint );          /* Debug Level */
432 
433 	}
434 	else if ( ! strcmp( optag, "dicom_coords" ) )
435 	{
436 	    str = PLUTO_get_string( plint );		/* Neg Extrema */
437 	    val = PLUTO_string_index(str, 2, grNY);
438 	    if ( val == 0 ) dicom_coords = 0;
439 	}
440 	else	/* illegal option? */
441 	{
442 	    sprintf( grMessage, "Error: pa_00\n"
443 		     "Unexpected option #%d: '%s'", opcnt, optag );
444 	    RETURN(grMessage);
445 	}
446 
447 	opcnt++;
448     }
449 
450     if ( ( out_rad > 0 ) && ( outfile == NULL ) )
451 	RETURN( "------------------------------------------------\n"
452                 "arguments : specify outfile to use output radius\n"
453                 "------------------------------------------------");
454 
455     if ( ! r_set_afni_s_from_dset( A, dset, debug ) )
456 	RETURN( "-------------------------------\n"
457                 "arguments : afni_s init failure\n"
458                 "-------------------------------");
459 
460     if ( ! init_maxima_s( M, A, outfile ) )
461 	RETURN("----------------------------------\n"
462                "MAXIMA_main: maxima_s init failure\n"
463                "----------------------------------");
464 
465     /* now fill any remaining parameters */
466     M->sval_style   = sval_style;
467     M->cutoff       = cutoff / A->factor[0];
468     M->min_dist     = min_dist;
469     M->out_rad      = out_rad;
470 
471     M->negatives    = negatives;
472     M->ngbr_style   = style;
473     M->quiet        = quiet;
474     M->true_max     = true_max;
475     M->dicom_coords = dicom_coords;
476     M->debug        = debug;
477 
478     gr_fac          = A->factor[0];
479 
480     if ( M->debug > 0 )
481     {
482 	if ( M->debug > 1 ) show_maxima_s( "plugin values applied ", M );
483 	fprintf(stderr,"  using sub-brick %d, factor %f (1/%f)\n",
484 		A->sub_brick, A->factor[0], 1/A->factor[0]);
485     }
486 
487     RETURN(NULL);
488 }
489 
490