00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*----------------------------------------------------------------------
00008  *
00009  *  plug_hemisub.c      - AFNI plugin to subtract hemispheres
00010  *
00011  *  $Log: plug_hemisub.c,v $
00012  *  Revision 1.6  2005/04/19 21:07:17  rwcox
00013  *  Cput
00014  *
00015  *  Revision 1.5  2004/01/07 19:50:37  rwcox
00016  *  Cput
00017  *
00018  *  Revision 1.4  2000/12/21 16:10:54  cox
00019  *  AFNI
00020  *
00021  *  Revision 1.2  2000/06/15 22:02:40  cox
00022  *  AFNI
00023  *
00024  *  Revision 1.1  1999/08/06 19:10:48  cox
00025  *  AFNI
00026  *
00027  * Revision 1.1  1998/09/04  19:59:54  rickr
00028  * Initial revision
00029  *
00030  *----------------------------------------------------------------------
00031  */
00032 
00033 #include "afni.h"
00034 
00035 typedef struct
00036 {
00037     int thresh_type;    /* none, pos_only, neg_only */
00038 } hemi_s;
00039 
00040 static char * process_data     ( THD_3dim_dataset *, hemi_s * );
00041 static char * process_as_floats( THD_3dim_dataset *, hemi_s * );
00042 
00043 #ifndef ALLOW_PLUGINS
00044 #  error "Plugins not properly set up -- see machdep.h"
00045 #endif
00046 
00047 /***********************************************************************
00048   plugin to find local extrema
00049 ************************************************************************/
00050 
00051 char * HEMISUB_main( PLUGIN_interface * );
00052 
00053 #define            NUM_T_OPTS   3
00054 
00055 static char      * thresh_opts[ NUM_T_OPTS ] =
00056                         { "any", "positives only", "negatives only" };
00057 static char        helpstring[] =
00058     "Hemisubtract - used to subtract one hemisphere from the other\n"
00059     ;
00060 
00061 
00062 
00063 /***********************************************************************
00064    Set up the interface to the user
00065 ************************************************************************/
00066 
00067 
00068 DEFINE_PLUGIN_PROTOTYPE
00069 
00070 PLUGIN_interface * PLUGIN_init( int ncall )
00071 {
00072     PLUGIN_interface * plint ;
00073 
00074     if( ncall > 0 ) return NULL ;  /* only one interface */
00075 
00076     /* create the new interface */
00077 
00078     plint = PLUTO_new_interface( "Hemi-subtract", "hemisphere subtraction",
00079                 helpstring, PLUGIN_CALL_VIA_MENU , HEMISUB_main );
00080 
00081     PLUTO_add_hint( plint,
00082         "from each voxel's value, subtract that of the reflected voxel" );
00083 
00084     PLUTO_set_sequence( plint , "z:Reynolds" ) ;
00085 
00086     /*-- first line of input: input dataset --*/
00087 
00088     PLUTO_add_option( plint, "Input" , "Input" , TRUE );
00089     PLUTO_add_hint( plint, "choose dataset for input" );
00090     PLUTO_add_dataset(plint, "Dataset" , ANAT_ALL_MASK , FUNC_ALL_MASK,
00091                                          DIMEN_3D_MASK | BRICK_SHORT_MASK );
00092 
00093     /*-- second line of input: prefix for output dataset --*/
00094 
00095     PLUTO_add_option( plint, "Output" , "prefix" , TRUE );
00096     PLUTO_add_hint( plint, "option: choose dataset prefix for output" );
00097     PLUTO_add_string( plint, "Prefix", 0, NULL, 19 );
00098 
00099     /*-- third line of input: threshold type option --*/
00100 
00101     PLUTO_add_option( plint, "Thresh Type", "Thresh Type", FALSE );
00102     PLUTO_add_string( plint, "Type", NUM_T_OPTS, thresh_opts, 0 );
00103 
00104     return plint;
00105 }
00106 
00107 
00108 /*----------------------------------------------------------------------
00109 **
00110 **  Main routine for this plugin (will be called from AFNI).
00111 **
00112 **----------------------------------------------------------------------
00113 */
00114 char * HEMISUB_main( PLUGIN_interface * plint )
00115 {
00116     THD_3dim_dataset * dset, * new_dset;
00117     MCW_idcode       * idc;
00118     hemi_s             hs = { 0 };
00119     char             * new_prefix;
00120     char             * ret_string = NULL;
00121     char             * tag;
00122 
00123 
00124     if ( plint == NULL )
00125         return  "------------------------\n"
00126                 "HEMISUB_main: NULL input\n"
00127                 "------------------------\n";
00128 
00129     PLUTO_next_option( plint );
00130     idc  = PLUTO_get_idcode( plint );
00131     dset = PLUTO_find_dset( idc );
00132 
00133     if( dset == NULL )
00134         return "-------------------------------\n"
00135                "HEMISUB_main: bad input dataset\n"
00136                "-------------------------------";
00137 
00138     DSET_load( dset );
00139 
00140     PLUTO_next_option( plint );
00141     new_prefix = PLUTO_get_string( plint );
00142     if ( ! PLUTO_prefix_ok( new_prefix ) )
00143         return  "------------------------\n"
00144                 "HEMISUB_main: bad prefix\n"
00145                 "------------------------\n";
00146 
00147     if ( ( new_dset = PLUTO_copy_dset( dset, new_prefix ) ) == NULL )
00148         return  "------------------------------------------\n"
00149                 "HEMISUB_main: failed to copy input dataset\n"
00150                 "------------------------------------------\n";
00151 
00152     tag = PLUTO_get_optiontag( plint );
00153     if ( tag && ! strcmp( tag, "Thresh Type" ) )
00154     {
00155         tag = PLUTO_get_string( plint );
00156         if ( tag != NULL )
00157             hs.thresh_type = PLUTO_string_index( tag, NUM_T_OPTS, thresh_opts );
00158     }
00159 
00160     if ( ret_string = process_data( new_dset, &hs ) )
00161         return  ret_string;
00162 
00163     if ( PLUTO_add_dset( plint, new_dset, DSET_ACTION_MAKE_CURRENT ) )
00164     {
00165         THD_delete_3dim_dataset( new_dset, False );
00166 
00167         return  "---------------------------------------\n"
00168                 "HEMISUB_main: failed to add new dataset\n"
00169                 "---------------------------------------\n";
00170     }
00171 
00172     return NULL;
00173 }
00174 
00175 
00176 /*----------------------------------------------------------------------
00177 **
00178 **  Subtract hemispheres.
00179 **
00180 **  Check if we need to create or change the factor.
00181 **
00182 **----------------------------------------------------------------------
00183 */
00184 static char *
00185 process_data( THD_3dim_dataset * dset, hemi_s * hs )
00186 {
00187     int     count, nx, ny, nz, cx, cx2;
00188     int     type, diff, floats = ( DSET_BRICK_FACTOR( dset, 0 ) != 0.0 );
00189     short * data, * sp, * sp2;
00190 
00191     nx = dset->daxes->nxx;
00192     ny = dset->daxes->nyy;
00193     nz = dset->daxes->nzz;
00194 
00195     type = hs->thresh_type;
00196 
00197     data = (short *)DSET_ARRAY( dset, 0 );
00198     for ( count = 0; ! floats && count < ny*nz; count++ )
00199     {
00200         sp  = data;
00201         sp2 = data + nx - 1;
00202 
00203         for ( cx = 0; cx < (nx+1)/2; cx++ )
00204         {
00205             if ( type == 1 )            /* positives only */
00206             {
00207                 if ( *sp < 0 )
00208                     *sp = 0;
00209                 if ( *sp2 < 0 )
00210                     *sp2 = 0;
00211             }
00212             else if ( type == 2 )       /* negatives only */
00213             {
00214                 if ( *sp > 0 )
00215                     *sp = 0;
00216                 if ( *sp2 > 0 )
00217                     *sp2 = 0;
00218             }
00219 
00220             diff = *sp - *sp2;
00221                                                   /* if out of short range */
00222             if ( ( diff > 32767 ) || ( diff < -32768 ) )
00223                 floats = 1;
00224             else
00225             {
00226                 *sp  = diff;
00227                 *sp2 = -diff;
00228             }
00229 
00230             sp++;
00231             sp2--;
00232         }
00233 
00234         data += nx;
00235     }
00236 
00237     if ( floats )
00238         return process_as_floats( dset, hs );
00239 
00240     return NULL;        /* success */
00241 }
00242 
00243 
00244 /*----------------------------------------------------------------------
00245 **
00246 **  Subtract hemispheres assuming we need floats.
00247 **
00248 **----------------------------------------------------------------------
00249 */
00250 static char *
00251 process_as_floats( THD_3dim_dataset * dset, hemi_s * hs )
00252 {
00253     int     count, cx, type = hs->thresh_type;
00254     int     nx, ny, nz, nvox;
00255     short * sp, * sdata;
00256     float * fdata, * fp, * fp2;
00257     float   factor, maxabs;
00258 
00259     nx   = dset->daxes->nxx;
00260     ny   = dset->daxes->nyy;
00261     nz   = dset->daxes->nzz;
00262     nvox = nx * ny * nz;
00263 
00264     sdata = (short *)DSET_ARRAY( dset, 0 );
00265 
00266     factor = DSET_BRICK_FACTOR( dset, 0 );
00267     factor = factor == 0.0 ? 1.0 : factor;
00268 
00269     /* first get the data into a float array */
00270 
00271     if ( ( fdata = (float *)malloc( nvox * sizeof( float ) ) ) == NULL )
00272         return  "------------------------------\n"
00273                 "paf: failed allocation of floats"
00274                 "------------------------------\n";
00275 
00276     fp = fdata;
00277     sp = sdata;
00278     for ( count = 0; count < nvox; count++ )
00279     {
00280         *fp = *sdata * factor;
00281 
00282         if ( ( type == 1 ) && ( *fp < 0 ) )
00283             *fp = 0;
00284         else if ( ( type == 2 ) && ( *fp > 0 ) )
00285             *fp = 0;
00286 
00287         fp++;
00288         sp++;
00289     }
00290 
00291     /* now make the subtraction as floats */
00292 
00293     for ( count = 0; count < ny*nz; count++ )
00294     {
00295         fp  = fdata + count * nx;
00296         fp2 = fp + nx - 1;
00297 
00298         for ( cx = 0; cx < (nx+1)/2; cx++ )
00299         {
00300             *fp  = *fp - *fp2;
00301             *fp2 = -*fp;
00302 
00303             fp++;
00304             fp2--;
00305         }
00306     }
00307 
00308     /* now make a new factor */
00309 
00310     maxabs = MCW_vol_amax( nvox, 1, 1, MRI_float, fdata );
00311 
00312     /* result is all zero, let the user worry */
00313     if ( maxabs != 0.0 )
00314     {
00315         factor = MRI_TYPE_maxval[MRI_short] /maxabs;        /* 32767? / maxabs */
00316     
00317         EDIT_coerce_scale_type( nvox, factor, MRI_float, fdata, MRI_short, sdata );
00318     
00319         DSET_BRICK_FACTOR( dset, 0 ) = factor == 0.0 ? 0.0 : 1.0 / factor;
00320     
00321         THD_load_statistics( dset );
00322     }
00323     free(fdata);
00324     return NULL;        /* success */
00325 }
00326