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 }

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 }

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 );         /* replace old point list with new */
00639     M->P = newP;
00640 
00641     return 1;
00642 }

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 };  /* for allocation speed */
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                                                 /* maybe a switch later */
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 }

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 }

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             /* do dicom coordinates from ijk, if requested */
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 }

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 ];               /* for speed */
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;          /* skip first plane */
01023     destp   = M->result + nxy;
01024 
01025     for ( cz = 0; cz < M->nz-2; cz++ )  /* and skip last plane (1->2) [rickr] */
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 }

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 }

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 }

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;           /* output will be 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 }

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 }

00202 {
00203    PLUGIN_interface * plint ;
00204 
00205    if( ncall > 0 ) return NULL ;  /* only one interface */
00206 
00207    /* create the new interface */
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    /*-- first line of input: input dataset --*/
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 ) ; /* new [rickr] */
00223 
00224 
00225    /*-- second line of input: prefix for output dataset --*/
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    /*-- third line of input: cutoff option --*/
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    /*-- fourth line of input: min_dist option --*/
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    /*-- fifth line of input: out_rad option --*/
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    /*-- sixth line of input: style option --*/
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    /*-- seventh line of input: negatives and true max options --*/
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    /*-- eighth line of input: true_max option --*/
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    /*-- ninth line of input: dicom_coords option --*/
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 }

01378 {
01379     short v1 = *( gr_orig_data + *(int *)p1 );
01380     short v2 = *( gr_orig_data + *(int *)p2 );
01381 
01382     /* more negative is greater AND reverse the order */
01383     if ( v1 > v2 )
01384         return 1;
01385     else if ( v1 < v2 )
01386         return -1;
01387 
01388 
01389     return 0;
01390 }

01357 {
01358     short v1 = *( gr_orig_data + *(int *)p1 );
01359     short v2 = *( gr_orig_data + *(int *)p2 );
01360 
01361     if ( v1 < v2 )              /* reverse order to get large numbers first */
01362         return 1;
01363     else if ( v1 > v2 )
01364         return -1;
01365 
01366     return 0;
01367 }

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     /* get AFNI inputs */
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 );        /* 2004 Feb 20 [rickr] */
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 );            /* Neg Extrema */
00404             val = PLUTO_string_index(str, 2, grNY);
00405             if ( val > 0 ) negatives = 1;
00406             str = PLUTO_get_string( plint );            /* True Max    */
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 );            /* No Text Out */
00430             val = PLUTO_string_index(str, 2, grNY);
00431             if ( val > 0 ) quiet = 1;
00432             debug = PLUTO_get_number( plint );          /* Debug Level */
00433             
00434         }
00435         else if ( ! strcmp( optag, "dicom_coords" ) )
00436         {
00437             str = PLUTO_get_string( plint );            /* Neg Extrema */
00438             val = PLUTO_string_index(str, 2, grNY);
00439             if ( val == 0 ) dicom_coords = 0;
00440         }
00441         else    /* illegal option? */
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     /* now fill any remaining parameters */
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 }

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;            /* global needed for sorting */
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 }

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 }

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;                 /* rickr - use sub-brick */
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 ];   /* just for speed */
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 /*    A->num_dsets++;   not using more than one */
01461 
01462     RETURN(1);
01463 }

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 }

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 }

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;         /* we're done */
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 }

00757 {
00758     double  total_x = 0.0, total_y = 0.0, total_z = 0.0;  /* weight*position */
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 )         /* no weighting necessary */
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 );      /* ~rounded average */
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 }

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     /* actually write a new dataset */
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 }