/*****************************************************************************
Major portions of this software are copyrighted by the Medical College
of Wisconsin, 1994-2000, and are released under the Gnu General Public
License, Version 2. See the file README.Copyright for details.
******************************************************************************/
/*----------------------------------------------------------------------
*
* plug_hemisub.c - AFNI plugin to subtract hemispheres
*
* $Log: plug_hemisub.c,v $
* Revision 1.6 2005/04/19 21:07:17 rwcox
* Cput
*
* Revision 1.5 2004/01/07 19:50:37 rwcox
* Cput
*
* Revision 1.4 2000/12/21 16:10:54 cox
* AFNI
*
* Revision 1.2 2000/06/15 22:02:40 cox
* AFNI
*
* Revision 1.1 1999/08/06 19:10:48 cox
* AFNI
*
* Revision 1.1 1998/09/04 19:59:54 rickr
* Initial revision
*
*----------------------------------------------------------------------
*/
#include "afni.h"
typedef struct
{
int thresh_type; /* none, pos_only, neg_only */
} hemi_s;
static char * process_data ( THD_3dim_dataset *, hemi_s * );
static char * process_as_floats( THD_3dim_dataset *, hemi_s * );
#ifndef ALLOW_PLUGINS
# error "Plugins not properly set up -- see machdep.h"
#endif
/***********************************************************************
plugin to find local extrema
************************************************************************/
char * HEMISUB_main( PLUGIN_interface * );
#define NUM_T_OPTS 3
static char * thresh_opts[ NUM_T_OPTS ] =
{ "any", "positives only", "negatives only" };
static char helpstring[] =
"Hemisubtract - used to subtract one hemisphere from the other\n"
;
/***********************************************************************
Set up the interface to the user
************************************************************************/
DEFINE_PLUGIN_PROTOTYPE
PLUGIN_interface * PLUGIN_init( int ncall )
{
PLUGIN_interface * plint ;
if( ncall > 0 ) return NULL ; /* only one interface */
/* create the new interface */
plint = PLUTO_new_interface( "Hemi-subtract", "hemisphere subtraction",
helpstring, PLUGIN_CALL_VIA_MENU , HEMISUB_main );
PLUTO_add_hint( plint,
"from each voxel's value, subtract that of the reflected voxel" );
PLUTO_set_sequence( plint , "z:Reynolds" ) ;
/*-- first line of input: input dataset --*/
PLUTO_add_option( plint, "Input" , "Input" , TRUE );
PLUTO_add_hint( plint, "choose dataset for input" );
PLUTO_add_dataset(plint, "Dataset" , ANAT_ALL_MASK , FUNC_ALL_MASK,
DIMEN_3D_MASK | BRICK_SHORT_MASK );
/*-- second line of input: prefix for output dataset --*/
PLUTO_add_option( plint, "Output" , "prefix" , TRUE );
PLUTO_add_hint( plint, "option: choose dataset prefix for output" );
PLUTO_add_string( plint, "Prefix", 0, NULL, 19 );
/*-- third line of input: threshold type option --*/
PLUTO_add_option( plint, "Thresh Type", "Thresh Type", FALSE );
PLUTO_add_string( plint, "Type", NUM_T_OPTS, thresh_opts, 0 );
return plint;
}
/*----------------------------------------------------------------------
**
** Main routine for this plugin (will be called from AFNI).
**
**----------------------------------------------------------------------
*/
char * HEMISUB_main( PLUGIN_interface * plint )
{
THD_3dim_dataset * dset, * new_dset;
MCW_idcode * idc;
hemi_s hs = { 0 };
char * new_prefix;
char * ret_string = NULL;
char * tag;
if ( plint == NULL )
return "------------------------\n"
"HEMISUB_main: NULL input\n"
"------------------------\n";
PLUTO_next_option( plint );
idc = PLUTO_get_idcode( plint );
dset = PLUTO_find_dset( idc );
if( dset == NULL )
return "-------------------------------\n"
"HEMISUB_main: bad input dataset\n"
"-------------------------------";
DSET_load( dset );
PLUTO_next_option( plint );
new_prefix = PLUTO_get_string( plint );
if ( ! PLUTO_prefix_ok( new_prefix ) )
return "------------------------\n"
"HEMISUB_main: bad prefix\n"
"------------------------\n";
if ( ( new_dset = PLUTO_copy_dset( dset, new_prefix ) ) == NULL )
return "------------------------------------------\n"
"HEMISUB_main: failed to copy input dataset\n"
"------------------------------------------\n";
tag = PLUTO_get_optiontag( plint );
if ( tag && ! strcmp( tag, "Thresh Type" ) )
{
tag = PLUTO_get_string( plint );
if ( tag != NULL )
hs.thresh_type = PLUTO_string_index( tag, NUM_T_OPTS, thresh_opts );
}
if ( ret_string = process_data( new_dset, &hs ) )
return ret_string;
if ( PLUTO_add_dset( plint, new_dset, DSET_ACTION_MAKE_CURRENT ) )
{
THD_delete_3dim_dataset( new_dset, False );
return "---------------------------------------\n"
"HEMISUB_main: failed to add new dataset\n"
"---------------------------------------\n";
}
return NULL;
}
/*----------------------------------------------------------------------
**
** Subtract hemispheres.
**
** Check if we need to create or change the factor.
**
**----------------------------------------------------------------------
*/
static char *
process_data( THD_3dim_dataset * dset, hemi_s * hs )
{
int count, nx, ny, nz, cx, cx2;
int type, diff, floats = ( DSET_BRICK_FACTOR( dset, 0 ) != 0.0 );
short * data, * sp, * sp2;
nx = dset->daxes->nxx;
ny = dset->daxes->nyy;
nz = dset->daxes->nzz;
type = hs->thresh_type;
data = (short *)DSET_ARRAY( dset, 0 );
for ( count = 0; ! floats && count < ny*nz; count++ )
{
sp = data;
sp2 = data + nx - 1;
for ( cx = 0; cx < (nx+1)/2; cx++ )
{
if ( type == 1 ) /* positives only */
{
if ( *sp < 0 )
*sp = 0;
if ( *sp2 < 0 )
*sp2 = 0;
}
else if ( type == 2 ) /* negatives only */
{
if ( *sp > 0 )
*sp = 0;
if ( *sp2 > 0 )
*sp2 = 0;
}
diff = *sp - *sp2;
/* if out of short range */
if ( ( diff > 32767 ) || ( diff < -32768 ) )
floats = 1;
else
{
*sp = diff;
*sp2 = -diff;
}
sp++;
sp2--;
}
data += nx;
}
if ( floats )
return process_as_floats( dset, hs );
return NULL; /* success */
}
/*----------------------------------------------------------------------
**
** Subtract hemispheres assuming we need floats.
**
**----------------------------------------------------------------------
*/
static char *
process_as_floats( THD_3dim_dataset * dset, hemi_s * hs )
{
int count, cx, type = hs->thresh_type;
int nx, ny, nz, nvox;
short * sp, * sdata;
float * fdata, * fp, * fp2;
float factor, maxabs;
nx = dset->daxes->nxx;
ny = dset->daxes->nyy;
nz = dset->daxes->nzz;
nvox = nx * ny * nz;
sdata = (short *)DSET_ARRAY( dset, 0 );
factor = DSET_BRICK_FACTOR( dset, 0 );
factor = factor == 0.0 ? 1.0 : factor;
/* first get the data into a float array */
if ( ( fdata = (float *)malloc( nvox * sizeof( float ) ) ) == NULL )
return "------------------------------\n"
"paf: failed allocation of floats"
"------------------------------\n";
fp = fdata;
sp = sdata;
for ( count = 0; count < nvox; count++ )
{
*fp = *sdata * factor;
if ( ( type == 1 ) && ( *fp < 0 ) )
*fp = 0;
else if ( ( type == 2 ) && ( *fp > 0 ) )
*fp = 0;
fp++;
sp++;
}
/* now make the subtraction as floats */
for ( count = 0; count < ny*nz; count++ )
{
fp = fdata + count * nx;
fp2 = fp + nx - 1;
for ( cx = 0; cx < (nx+1)/2; cx++ )
{
*fp = *fp - *fp2;
*fp2 = -*fp;
fp++;
fp2--;
}
}
/* now make a new factor */
maxabs = MCW_vol_amax( nvox, 1, 1, MRI_float, fdata );
/* result is all zero, let the user worry */
if ( maxabs != 0.0 )
{
factor = MRI_TYPE_maxval[MRI_short] /maxabs; /* 32767? / maxabs */
EDIT_coerce_scale_type( nvox, factor, MRI_float, fdata, MRI_short, sdata );
DSET_BRICK_FACTOR( dset, 0 ) = factor == 0.0 ? 0.0 : 1.0 / factor;
THD_load_statistics( dset );
}
free(fdata);
return NULL; /* success */
}
syntax highlighted by Code2HTML, v. 0.9.1