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 #include "afni.h"
00008 
00009 #ifndef ALLOW_PLUGINS
00010 #  error "Plugins not properly set up -- see machdep.h"
00011 #endif
00012 
00013 /***********************************************************************
00014   Simple plugin to compute dot product of two datasets.
00015 ************************************************************************/
00016 
00017 static char helpstring[] =
00018   " Purpose: Compute correlation of two 3D datasets\n"
00019   " Inputs:\n"
00020   " Dataset 1   = 3D dataset that must already be in memory\n"
00021   " Dataset 2   = 3D dataset that must already be in memory\n"
00022   " Remove Mean = If this option is toggled, then the means\n"
00023   "                of each dataset will be removed before\n"
00024   "                the correlations are computed.\n"
00025   " The output is popped up into a window on the screen.\n"
00026   "Author -- RW Cox"
00027 ;
00028 
00029 /*---------- prototypes for internal routines ----------*/
00030 
00031 static char * DOT_main( PLUGIN_interface * ) ;
00032 
00033 static double DSET_cor( int , THD_3dim_dataset * , THD_3dim_dataset * ) ;
00034 
00035 /***********************************************************************
00036    Set up the interface to the user
00037 ************************************************************************/
00038 
00039 DEFINE_PLUGIN_PROTOTYPE
00040 
00041 PLUGIN_interface * PLUGIN_init( int ncall )
00042 {
00043    PLUGIN_interface * plint ;
00044 
00045    if( ncall > 0 ) return NULL ;  /* only one interface */
00046 
00047    /*-- set titles and call point --*/
00048 
00049    plint = PLUTO_new_interface( "3D Correlation" , "3D Dataset Correlation" , helpstring ,
00050                                  PLUGIN_CALL_VIA_MENU , DOT_main  ) ;
00051 
00052    PLUTO_set_sequence( plint , "A:afniinfo:dset" ) ;
00053 
00054    PLUTO_add_hint( plint , "3D Dataset Correlation" ) ;
00055 
00056    /*-- first line of input: Dataset --*/
00057 
00058    PLUTO_add_option( plint , "Dataset" , "Dataset" , TRUE ) ;
00059    PLUTO_add_dataset(plint , "# 1" ,
00060                                     ANAT_ALL_MASK , FUNC_ALL_MASK ,
00061                                     SESSION_ALL_MASK |
00062                                     DIMEN_3D_MASK    | BRICK_ALLREAL_MASK ) ;
00063 
00064    /*-- second line of input: Dataset --*/
00065 
00066    PLUTO_add_option( plint , "Dataset" , "Dataset" , TRUE ) ;
00067    PLUTO_add_dataset(plint , "# 2" ,
00068                                     ANAT_ALL_MASK , FUNC_ALL_MASK ,
00069                                     DIMEN_3D_MASK | BRICK_ALLREAL_MASK ) ;
00070 
00071    /*-- third line of input: Remove Mean option --*/
00072 
00073    PLUTO_add_option( plint , "Remove Mean" , "Remove Mean" , False ) ;
00074 
00075    return plint ;
00076 }
00077 
00078 /***************************************************************************
00079   Main routine for this plugin (will be called from AFNI).
00080 ****************************************************************************/
00081 
00082 static char * DOT_main( PLUGIN_interface * plint )
00083 {
00084    MCW_idcode * idc ;
00085    THD_3dim_dataset * xset , * yset ;
00086    char * tag ;
00087    int demean ;
00088    double dcor ;
00089    char str[256] ;
00090 
00091    /*--------------------------------------------------------------------*/
00092    /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
00093 
00094    PLUTO_next_option(plint) ;                             /* go to next input line */ 
00095    idc  = PLUTO_get_idcode(plint) ; /* get 1st dataset item */
00096    xset = PLUTO_find_dset(idc) ;                   /* get ptr to dataset */
00097    if( xset == NULL )
00098       return "**********************\n"
00099              "Cannot find Dataset #1\n"
00100              "**********************"  ;
00101 
00102    PLUTO_next_option(plint) ;
00103    idc  = PLUTO_get_idcode(plint) ;
00104    yset = PLUTO_find_dset(idc) ;
00105    if( yset == NULL )
00106       return "**********************\n"
00107              "Cannot find Dataset #2\n"
00108              "**********************"  ;
00109 
00110    /*-- The first two input lines (processed above) were mandatory.
00111         The next one is optional.  Check to see if it is present,
00112         and if it has the right name.
00113         If so, then mark its presence in the variable "demean".   --*/
00114 
00115    tag    = PLUTO_get_optiontag(plint) ;
00116    demean = ( tag != NULL && strcmp(tag,"Remove Mean") == 0 ) ;
00117 
00118    /*------------------------------------------------------*/
00119    /*---------- At this point, the inputs are OK ----------*/
00120 
00121    /*-- do the actual work --*/
00122 
00123    dcor = DSET_cor( demean , xset , yset ) ;
00124 
00125    if( dcor < -1.0 )
00126       return "*********************************\n"
00127              "Error while computing correlation\n"
00128              "*********************************"  ;
00129 
00130    /*-- put the output to the screen --*/
00131 
00132    sprintf(str , "    Dataset %s\n"
00133                  "and Dataset %s\n\n"
00134                  "Correlation = %g\n"
00135                  "%s" ,
00136            DSET_FILECODE(xset) , DSET_FILECODE(yset) ,
00137            dcor ,
00138            (demean) ? "[mean removed]" : " " ) ;
00139 
00140    PLUTO_popup_message( plint , str ) ;
00141 
00142    return NULL ;  /* null string returned means all was OK */
00143 }
00144 
00145 /*----------------------------------------------------------------------------
00146   Compute the correlation between two datasets.  Returns a number
00147   less than -1.0 if an error occurs.
00148 ------------------------------------------------------------------------------*/
00149 
00150 static double DSET_cor( int demean, THD_3dim_dataset * xset, THD_3dim_dataset * yset )
00151 {
00152    double sumxx , sumyy , sumxy , tx,ty , dxy , xbar,ybar ;
00153    void  *  xar , *  yar ;
00154    float * fxar , * fyar ;
00155    int ii , nxyz , ivx,ivy , itypx,itypy , fxar_new,fyar_new ;
00156 
00157    /*-- check datasets for conformity in dimensions --*/
00158 
00159    nxyz = xset->daxes->nxx * xset->daxes->nyy * xset->daxes->nzz ;
00160 
00161    if(  yset->daxes->nxx * yset->daxes->nyy * yset->daxes->nzz != nxyz )
00162       return -666.0 ;
00163 
00164    /*-- load the first dataset into memory --*/
00165 
00166    DSET_load( xset ) ;
00167 
00168    ivx   = DSET_PRINCIPAL_VALUE(xset) ;  /* most important place */
00169    itypx = DSET_BRICK_TYPE(xset,ivx) ;   /* type of data stored here */
00170    xar   = DSET_ARRAY(xset,ivx) ;        /* get the array */
00171    if( xar == NULL ){
00172       DSET_unload(xset) ;  /* free memory if an error happened */
00173       return -666.0 ;
00174    }
00175 
00176    /*-- if the array is not floating point, make a floating point copy --*/
00177 
00178    if( itypx == MRI_float ){
00179       fxar = (float *) xar ; fxar_new = 0 ;
00180    } else {
00181       fxar = (float *) malloc( sizeof(float) * nxyz ) ; fxar_new = 1 ;
00182       EDIT_coerce_type( nxyz , itypx,xar , MRI_float,fxar ) ;
00183       DSET_unload( xset ) ;  /* don't need this in memory anymore */
00184    }
00185 
00186    /*-- do the same for the second dataset --*/
00187 
00188    DSET_load( yset ) ;
00189 
00190    ivy   = DSET_PRINCIPAL_VALUE(yset) ;  /* most important place */
00191    itypy = DSET_BRICK_TYPE(yset,ivy) ;   /* type of data stored here */
00192    yar   = DSET_ARRAY(yset,ivy) ;        /* get the array */
00193    if( yar == NULL ){
00194       if( fxar_new ) free(fxar) ; else DSET_unload(xset) ;  /* free memory */
00195       DSET_unload(yset) ;                                   /* if an error */
00196       return -666.0 ;
00197    }
00198 
00199    /*-- if the array is not floating point, make a floating point copy --*/
00200 
00201    if( itypy == MRI_float ){
00202       fyar = (float *) yar ; fyar_new = 0 ;
00203    } else {
00204       fyar = (float *) malloc( sizeof(float) * nxyz ) ; fyar_new = 1 ;
00205       EDIT_coerce_type( nxyz , itypy,yar , MRI_float,fyar ) ;
00206       DSET_unload( yset ) ;  /* don't need this in memory anymore */
00207    }
00208 
00209    /*-- if needed, compute the mean of each dataset --*/
00210 
00211    xbar = ybar = 0.0 ;
00212    if( demean ){
00213       for( ii=0 ; ii < nxyz ; ii++ ){
00214          xbar += fxar[ii] ;
00215          ybar += fyar[ii] ;
00216       }
00217       xbar /= nxyz ;
00218       ybar /= nxyz ;
00219    }
00220 
00221    /*-- now compute the quadratic sums --*/
00222 
00223    sumxx = sumyy = sumxy = 0.0 ;
00224    for( ii=0 ; ii < nxyz ; ii++ ){
00225       tx = (fxar[ii]-xbar) ; ty = (fyar[ii]-ybar) ;
00226       sumxx += tx * tx ;
00227       sumyy += ty * ty ;
00228       sumxy += tx * ty ;
00229    }
00230 
00231    /*-- free up arrays --*/
00232 
00233    if( fxar_new ) free(fxar) ; else DSET_unload(xset) ;
00234    if( fyar_new ) free(fyar) ; else DSET_unload(yset) ;
00235 
00236    /*-- compute correlation --*/
00237 
00238    dxy = sumxx * sumyy ; if( dxy <= 0.0 ) return -666.0 ;
00239    dxy = sumxy / sqrt(dxy) ;
00240    return dxy ;
00241 }